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SECTION  1.  INTRODUCTION 


This  report  documents  the  research  and  development  work  on  a  computer 
simulation  of  a  pyroelectric  thermal  focal  plane  array  by  Nichols  Research 
Corporation.  This  model  is  completely  operational  and  described  in  software 
documents  [1,2].  The  present  description  is  intended  primarily  for  individuals 
who  would  like  details  of  the  model's  formulation,  range  of  applicability, 
assumptions,  approximations  and  results.  These  topics  are  summarized  in 
Section  2. 

Section  3  of  this  document  covers  the  results  of  the  literature  survey. 
Main  outputs  of  the  survey  are  an  annotated  bibliography  which  is  published  in 
a  separate  volume  [3],  and  tables  of  pyroelectric  system  parameters  given  in 
Appendices  A  and  B. 

Section  4  presents  the  details  of  the  formulation  of  thermal  and  electrical 
parts  of  the  model.  Included  in  the  thermal  formulation  are  the  conductive, 
radiative  and  convective  exchange  within  and  between  the  detectors.  The 
electrical  formulation  includes  relating  detector  voltages  to  their  temperature 
changes;  the  various  transfer  functions  involved  in  signal  readout,  amplification 
and  processing;  and  the  formulation  of  system  performance  parameters. 

The  numerical  procedures  that  were  used  in  the  computer  code  to  evaluate 
its  outputs  are  described  in  Section  5.  Procedures  for  the  preprocessor, 
processor  and  postprocessor  are  covered  in  Sections  5.1,  5.2  and  5.3  respectively. 

A  wide  variety  of  validations  of  the  code  were  run  against  exact  analytical 
solutions  to  theoretical  situations  that  test  various  parts  of  the  code  (see  Section 
6).  Included  in  these  tests  are  different  boundary  conditions  and  configurations 
to  verify  the  conductive,  radiative  and  convective  procedures  used.  In  Section 
7,  we  evaluated  the  responsivity  and  crosstalk  parameters  of  typical  systems 
irradiated  by  a  temporally  chopped  line  source,  and  compared  the  results  to 
similar  analyses  on  a  typical  state-of-the-art  system.  Finally,  the  conclusions 
and  recommendations  of  the  study  are  discussed  in  Sections  8  and  9, 
respectively. 


SECTION  2.  SUMMARY 


The  development  of  the  computer  simulation  for  pyroelectric  thermal  focal 
plane  arrays  was  based  on  a  thorough  survey  of  the  available  literature  and  on 
discussions  with  several  system  developers.  The  model  can  handle  several 
general  system  configurations  with  arbitrary  thermal  and  electrical  properties. 
It  treats  the  transient  three-dimensional  thermal  response  of  the  detector 
material  and  the  associated  electrodes,  insulation  and  interconnectors,  as  they 
are  influenced  by  arbitrarily  specified,  chopped  thermal  radiation  and  other 
environmental  heat  transfer  mechanisms  (see  Figure  1). 

A  finite  difference  numerical  solution  was  developed  to  consider  general 
conductive,  radiative  and  convective  heat  transfer  within  and  between  pyro¬ 
electric  detector  elements.  The  user  can  specify  nonreticulated  or  fully 
reticulated  focal  planes.  The  fully  reticulated  gaps  can  be  backfilled  with 
arbitrary  insulation  or  gas,  or  completely  evacuated.  Arbitrary  focal  plane 
configurations  can  be  treated  with  arbitrary  internal  thermal  node  resolution 
within  practical  computer  run  time  and  storage  limitations.  The  signal 
production,  readout,  and  processing  mechanisms  and  their  associated  noise 
sources  are  fully  modeled  with  currently  accepted  theoretical  expressions. 

Various  system  performance  parameters  such  as  the  modulation  transfer 
function  (MTF),  detectivity  (D*),  responsivity,  and  noise  equivalent  power  are 
automatically  evaluated  for  the  simulated  sensor.  The  frequency  response  of 
several  sensors  modeled  by  Honeywell  Electro-Optics  Division  (HEOD)  are  in 
general  agreement,  but  about  a  factor  of  two  below  those  predicted  by  the 
present  model. 

A  variety  of  validations  to  the  thermal  model  were  completed  using 
different  exact  solutions  such  as  one-dimensional  heat  transfer  with  a  number  of 
boundary  conditions.  These  computer  runs  generally  checked  to  within  one 
percent  of  exact  solutions  for  pertinent  situations.  Comparisons  with  modeled 
crosstalk  results  published  by  HEOD  show  agreement  to  within  about  5  percent. 


The  code  was  designed  to  be  user-friendly  and  highly  modular  to  facilitate 
modifications  and  upgrades.  It  is  thoroughly  documented  both  internally  (by 
comment  cards  and  the  use  of  meaningful  names  for  code  parameters,  functions, 
and  subroutines),  and  externally  (by  a  User's  Manual  [1]  and  a  Software 
Maintenance  Manual  [2]). 

The  major  development  efforts  on  the  TFPA  model  since  the  Interim 
Report  [4]  were  on  the  radiative  and  convective  heat  transfer  between  detector 
elements,  the  relationship  between  detector  temperature  voltage,  and  a 
postprocessor  that  treats  all  readout  and  signal  processing  operations  and  their 
related  noise  contributions.  It  was  found  that  radiative  and  convective  transfer 
generally  has  a  negligible  effect  on  MTF  and  crosstalk. 
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FIGURE  1.  THERMAL  NODE  STRUCTURE,  COORDINATE  SYSTEM  AND  NOTATION. 


SECTION  3.  LITERATURE  SURVEY 


The  objective  of  this  task  was  to  gather  all  available  documentation  on 
pyroelectric  detectors  with  special  concern  toward  material  properties,  and 
treatments  of  signal  and  noise  processes  within  the  detector.  Our  approach  was 
to  conduct  computerized  searches  at  the  "DIALOG"  information  storage  and 
retrieval  system  and  the  ERIM/IR1A  Center  document  library.  Other  major 
sources  were  the  Night  Vision  and  Electro-Optics  Laboratory  (NVEOL),  and 
references  from  collected  documents.  Of  the  original  collection  of  221 
documents  that  were  referenced  in  our  Interim  Report  [4],  a  group  of  150  of  the 
most  pertinent  was  annotated  and  published  in  a  separate  volume  [3].  In 
addition,  several  of  the  most  promising  documents  of  that  group  were  reviewed 
in  detail  to  determine  the  most  appropriate  thermal  and  electrical  modeling 
procedures.  These  are  referenced  as  they  are  discussed  later  in  this  report 
(especially  Section  4). 

All  of  the  original  list  of  documents  were  also  reviewed  for  information  on 
material  properties  of  pyroelectric  materials.  A  table  of  pyroelectric  para¬ 
meters  (pyroelectric  constant,  dielectric  constant,  and  dielectric  loss  tangent) 
and  thermal  properties  (thermal  conductivity,  specific  heat  per  unit  mass,  mass 
density  per  unit  volume,  and  thermal  diffusivity)  are  given  for  38  different 
materials  (see  Appendix  A).  (Not  all  the  properties  were  available  for  all  the 
materials.)  The  source  of  these  properties  are  fully  referenced.  A  com¬ 
plementary  table  of  thermal  properties  is  given  (see  Appendix  B)  for  16  different 
typical  insulation  and  separation  materials  used  in  a  thermal  focal  plane  array. 
These  sources  are  also  fully  referenced. 


SECTION  4.  FORMULATION 


The  objective  of  this  section  is  to  fully  document  NRC's  formulation  of  the 
thermal  and  electrical  equations  that  are  amenable  to  computer  solution.  The 
approximations  invoked  are  fully  described  and  the  limits  of  applicability  of  the 
solution  indicated.  Most  of  this  development  is  not  new  but  rather  a  compilation 
of  existing  treatments  which  will  be  fully  referenced.  Subsections  4.1  and  4.2 
cover  the  thermal  and  electrical  formulations,  respectively. 

4.1  THERMAL  FORMULATION 

The  volumetric  averaged  temperatures  of  the  individual  detector  elements 
of  the  array  are  required  as  a  function  of  time  in  order  to  predict  the  detector 
output  voltage.  The  model  must  be  capable  of  treating  rather  general 
geometrical  configurations  (see  Figure  1)  for  the  elements  with  arbitrary 
thermal  properties  under  arbitrary  thermal  boundary  conditions.  We  have 
retained  the  three-dimensional  capability  in  order  to  handle  general  detector 
arrays  under  any  irradiance  patterns. 

The  usual  three-dimensional  heat  conduction  equation  for  an  isotropic 
material  with  constant  thermal  properties  is: 


where  T  is  the  temperature  of  the  material  as  a  function  of  time  t  and  the  three 
orthogonal  spatial  coordinates  (x,y,z),  and  *  is  the  thermal  diffusivity  of  the 
material.  The  boundary  conditions  we  employ  are  (1)  a  uniform  initial 
temperature  distribution: 


and  (2)  the  radiative,  convective,  and  conductive  heating  or  cooling  specifica«.iow 
at  the  exterior  surfaces. 


This  second  set  of  boundary  conditions  depends  on  the  details  of  the  system 
configuration  being  modeled.  The  irradiated  surface  of  the  detector  (z=0)  has 
a  net  radiative  gain  per  unit  area  of: 


ae<x,y,t)  +  eo 


T(x,y,0, 


(3) 


where  E  is  the  spatial  irradiance  difference  from  a  blackbody  environment  at  a 
uniform  background  temperature  T&,  a  and  e  are  the  absorptance  and  emittance 
of  the  irradiated  surface,  and  a  is  the  Stefan- Boltzmann  constant. 

The  convective  heat  transfer  is  generally  given  by: 

*cnv  "  hITg  *  T(x,y,0,t))  (4) 

where  h  is  a  correlation  constant  and  the  ambient  gas  is  taken  to  have  a 
constant  uniform  temperature  Tg.  To  evaluate  the  constant  h,  we  used  published 
correlation  results  for  natural  convection  to  surfaces  in  air  [51.  These 
correlations  are  given  for  different  orientations  and  thermal  conditions  (such  as 
vertical  and  horizontal  plates  hotter  or  cooler  than  the  ambient  gas).  We  have 
used  the  correlation  for  vertical  plates  since  the  orientation  of  the  surfaces  is 
unknown,  but  vertical  orientation  is  most  likely. 

hL/kf  =  0.555  (GrPr)l/3  (5) 

where  L  is  the  characteristic  length  (taken  to  be  the  length  of  the  detector),  kf 
is  the  thermal  conductivity  of  the  gas  film  over  the  surface,  Gr  is  the  Grashoff 
number  and  Pr  is  the  Prandtl  number.  Typical  room  temperature  values  of  kf 
and  Pr  were  taken  from  Reference  5.  The  Grashoff  number  is  given  by: 


where  Bg/v2  is  given  by  Kreith  [5].  Using  a  nominal  characteristic  length  (L  = 
10“4  m),  we  find: 

h  =  1.71  aT1^  J'S-i-m"2*^1  (7) 


so 


*Qcnv  =  1*71  (Tg  -  T(x,y,0,t)]  4/3  (8) 

where  the  units  of  q  are  .T^rn"2,  when  T  is  in  K. 

The  boundary  condition  at  the  other  surfaces  of  the  detector  depend  on  the 
details  of  the  array  configuration.  For  example,  one  option  of  the  code  allows 
for  insulation  material  between  adjacent  detectors  in  the  array.  In  this  case,  we 
neglect  radiative  and  convective  heat  transfer  from  these  surfaces,  but  consider 
instead  the  conduction  of  heat  through  the  insulation.  The  boundary  conditions 
simply  match  the  energy  flow  from  one  surface  to  the  energy  flow  to  the 
adjacent  surface.  That  is,  the  heat  leaving  the  side  of  the  detector  must  equal 
the  heat  entering  the  insulator  at  the  same  point.  Other  boundary  conditions  are 
specified  similarly. 

In  order  to  solve  the  heat  conduction  equation  (see  Equation  1)  numerically 
on  the  computer,  we  expand  it  in  finite  difference  form  [6J: 

C(i,j,k)[T(i,j,k,t  +  At)  -  T(i,J,k.t)]  - 


1  J  K 

I  I  I  K(i,j,k; 
iT»l  JT-l  kT-l 


i'.j'.k’)  [T<i\j\k\t)  -  T(i,J ,k,t) ]  At 


(9) 


where  C  is  the  heat  capacity  of  the  node  and  the  K's  are  the  thermal 
conductances  for  the  paths  to  the  adjacent  nodes  (indicated  by  the  primes).  The 
left  side  represents  the  energy  increase  of  the  node  in  terms  of  its  temperature 


rise  during  the  time  period,  and  the  right  represents  the  heat  conducted  into  the 
node  based  on  temperatures  of  the  previous  time  t.  As  was  discussed  earlier  [4], 
we  selected  the  explicit  treatment  for  programming  expediency,  but  would 
recommend  investigation  of  the  implicit  procedures  such  as  Crank-Nicholson  and 
others,  that  assure  convergence. 

In  this  equation,  the  temperature  of  the  node  under  consideration  at  the 
previous  calculational  time  T(i,j,k,t),  has  a  coefficient  (1  -  t/x),  where  the 
characteristic  time  is: 


/[  “i  nj  nk 

T  =  C(i.j.k)/  l  l  l  K(i,j,k;i\j\k') 

/  Li'-l  j’-l  k'-l 


This  characteristic  time  provides  the  criterion  for  absolute  stability.  It  can  be 
seen  by  inspection  that  when  the  time  step  At  is  greater  than  this  parameter, 
the  dependency  on  the  temperature  of  the  node  at  the  previous  calculational 
time  becomes  negative,  implying  instability.  The  preprocessor  of  the  simulation 
effectively  finds  the  largest  allowable  value  of  x  (i,j,k)  for  all  directions  guiding 
the  user  to  select  a  proper  value  of  the  calulational  time  step  At  to  assure 
calculational  stability. 

The  top  surface  boundary  condition  involves  the  optical  properties  (of )  of 
the  absorbing  material.  We  are  assuming  here  that  the  focused  radiation  (E) 
would  be  filtered  so  that  the  absorptance  (a)  should  be  evaluated  for  the 
bandpass  of  the  instrument  and  for  the  directions  of  the  incident  radiation.  The 
emittance  (e)  should  be  hemispherical  and  integrated  over  all  wavelengths.  We 
have  assumed  in  our  calculations  to  date  that  a  and  e  are  equal  to  a  nominal 
value  given  in  Appendix  A,  although  arbitrary  values  can  also  be  assigned  to 
each.  At  present,  we  can  select  the  effective  background  temperature,  Tp, 
independent  of  the  ambient  gas  temperature  Tg  to  accommodate  more  general 
background-atmospheric-operating  conditions.  The  chopper  irradiance  can  be 
specified  in  a  general  manner  to  account  for  the  arbitrary  conditions. 


The  boundary  conditions  on  the  sides  of  the  detectors  were  also  made  as 
general  as  possible.  The  user  can  specify  a  solid  material  between  the  detector 
elements  having  any  real  value  of  thermal  conductivity  (k)  desired  independent 
of  the  materials  below  the  detectors.  Thus,  any  theoretical  (e.g„  a  perfect 
insulation  with  k=0),  or  realizable  conductivity  such  as  for  any  typical  insulators 
(see  Appendix  B),  can  be  handled.  Non-retieulated  detectors  or  fully  reticulated 
detectors  can  also  be  treated  by  proper  selection  of  material  properties. 
(Reticulated  detectors  are  separated  by  a  gap  to  reduce  heat  transfer  between 
adjacent  detector  elements.) 

Stagnant  gas,  such  as  may  be  found  in  the  spaces  between  detectors,  and 
between  detectors  and  the  substrate  material  can  also  be  modeled  with  a 
different  value  of  conductivity  and  heat  capacity.  It  is  difficult  to  model  gas 
in  this  way,  however,  since  extremely  small  time  steps  are  required  for  stability 
(see  Equation  10). 

Convective  and/or  radiative  boundary  conditions  can  be  specified  for  the 
sides  and  bottoms  of  the  detectors.  The  convective  heat  transfer  coefficient  for 
these  surfaces  has  default  values  based  on  the  same  correlations  as  the  ones  used 
for  the  top  surface.  This  can  be  set  to  any  value  (including  zero)  for  modeling 
evacuated  or  partially  evacuated  systems. 

The  radiative  transfer  between  detectors,  and  between  detectors  and  the 
substrate  is  treated  assuming  that  only  the  two  directly  opposing  nodes  are 
directly  controlling  the  heat  transfer.  That  is,  we  consider  parallel  rectangular 
surfaces  of  equal  size  (representing  the  two  opposing  nodes)  separated  by  a 
distance,  D.  All  other  surfaces  of  the  radiative  enclosure  are  assumed  to  absorb 
just  as  much  heat  as  they  emit  (so-called  "refractory  surfaces").  For  this 
situation,  the  net  heat  transfer  from  Surface  1  to  Surface  2  is  [7J: 

*12  ■  »1>  (T,‘  -  Tj*)  <11> 

where  the  upper  case  Q  implies  the  total  heat  flux  over  A,  the  area  of  either 
surface.  The  radiation  factor  is  given  for  diffuse,  gray  surfaces  as: 


*12  "  <1^«1  ♦  »/c2  «■  1/Fj,)-1 


(12) 


The  geometrical  radiation  exchange  factor  is  given  by  Kreith  for  equal 
parallel  rectangles  with  width  W  and  length  L,  and  separated  by  a  distance  D  [5]: 


r,,.-  ~  In 
12  luv 


*  u(l+v  ) 


-1  — vcl 

1 1  l-Ki  4v  l<l4t.  )*J 

ir2)**  tan  1  - Vt  "  u  tan  *(*»)  “  v  tan  *(*) 

[a*  )*] 


where  u  =  D/L  and  v  =  D/W.  The  radiative  exchange  between  the  bottom  of  the 
detector  and  the  substrate  is  evaluated  the  same  way  with  the  area  in  Equation 
11  taken  to  be  the  bottom  surface  area  of  the  detector  node.  We  have  found, 
however,  that  for  the  low  temperature  and  small  temperature  differences 
encountered  in  this  problem  that  the  radiative  exchange  is  negligible  compared 
to  convection  or  conduction. 

The  conductive  heat  transfer  boundary  condition  for  the  bottom  of  the 
detector  considers  a  pedestal  exchanging  heat  with  the  substrate.  This  pedestal 
can  be  surrounded  by  a  second  insulation  material  at  the  user's  option.  Thus, 
the  thermal  properties  (conductivity  and  heat  capacity)  of  the  pedestal  and  two 
different  insulation  materials  can  be  specified  separately. 


It  should  be  noted  that  our  lower  boundary  condition  assumes  a  constant 
temperature  (i.e.,  thermally  massive)  substrate.  Texas  Instruments  concurs  with 
these  approximations  for  their  system.  On  the  other  hand,  HEOD's  thermal 
model  [8]  considers  a  thermal  conduction  path  between  the  detectors  through  the 
substrate. 


4.2  ELECTRICAL  MODEL 


The  thermal  model  provides  a  means  whereby  the  temperature  of  the 
array,  as  a  function  of  position  and  time,  can  be  determined  as  a  result  of  the 
input  irradiance  pattern  upon  it.  These  temperature  changes  are  translated  into 
voltages  within  the  array,  and  passed  through  a  readout  mechanism  to  be 
ultimately  analyzed  or  displayed.  The  electrical  portion  of  the  model  must 
describe  the  temperature  change  to  voltage  conversion  process,  and  the  circuitry 
used  to  extract  voltages  from  the  array. 

In  addition,  electronic  noise,  generated  in  the  array  elements,  must  also  be 
modeled.  These  noises  are  amplified  together  with  the  noises  generated  within 
the  readout  structure  to  produce  variations  in  the  array  output  which  are 
uncharacteristic  of  the  innut  scene. 

Once  the  electrical  characteristics  and  the  noises  are  modeled,  the 
performance  of  the  array  system  can  be  modeled  in  terms  of  parameters 
traditionally  used  to  describe  the  performance  of  real  devices,  such  as  MTF, 
noise  equivalent  power  (NEP),  responsivity,  and  specific  detectivity  (D*).  The 
ultimate  usefulness  of  this  code  comes  from  its  ability  to  model  both  the 
thermal  and  electrical  characteristics  of  thermal  arrays.  In  this  way,  the 
expected  performance  of  arrays  can  be  assessed  in  terms  of  their  characteristics 
without  the  requirement  to  actually  build  and  test  an  example. 

General  Temperature  Change  to  Voltage  Conversion 

The  transduction  of  temperature  change  to  voltage  in  the  subject  arrays 
occurs  through  the  pyroelectric  effect.  The  spatial  temperature  variations  in 
the  array  are  due  to  spatial  variations  in  the  radiative  heat  load.  For  a  given 
1/2  chopping  cycle  time  Az,  a  given  detector  element  (M,N)  will  experience  an 
effective  temperature  change  AT.  A  surface  electric  charge  is  generated  on  the 
detector  electrodes  through  the  pyroelectric  effect  and  the  magnitude  of  the 
charge  change  is  given  by: 


where  : 


Aq  =  charge  change  in  coulombs 
p  =  pyroelectric  coefficient  in  coulombs/K/cm2 
AT  =  temperature  change  in  K 
d  =  detector  area  in  cm2 

Another  mode  of  operation  that  has  proved  successful  is  the  use  of  an 
electrically  biased  dielectric  that  undergoes  a  dielectric  constant  change  as  a 
function  of  temperature  (T.I.  system).  The  equations  that  pertain  to  this  system 
are: 


qi  =  CjV  at  temperature  Tj 
<\2  -  C2V  at  temperature  T2 


AQ  =  <Jl  "  d2  =  v<cl  '  c2>  =  Ko<q  -  ®2>  =  KGAe  (15) 


If  we  define  a  term  pQ  =  de/dT  then: 

q  =  KoPoAT  (16) 

and  the  analytic  equations  are  essentially  the  same  with  slightly  different 
physical  meanings  attributed  to  the  figures-of-merit.  The  array  does  not 
transform  the  irradiance  input  into  temperature  (or  voltage)  perfectly,  nor  does 
the  readout  device  transfer  the  voltages  out  of  the  array  perfectly,  so  that  an 
ideal  display  will  show  spatial  distortions  due  to  thermal  spread  between 
elements  and  to  electrical  transfer  inefficiencies.  In  addition,  electrical  noise 
generated  within  the  array  elements  and  in  the  input  circuitry  will  be  added  to 
the  array  element  voltages.  These  noises  are  responsible  for  ultimately  limiting 
the  ability  of  the  array  to  detect  spatial  variations  in  irradiance. 

The  ultimate  goal  of  the  electrical  modeling  is  to  characterize  the 
performance  of  the  array.  We  wish  to  model  its  spatial  resolving  capability,  by 
use  of  the  MTF,  its  noise  limited  sensitivity,  by  use  of  the  noise  equivalent 
irradiance  or  flux  density  (NEFD),  and  its  irradiance-to-voltage  transfer 
function,  its  responsivity. 


Performance  modeling  of  the  array  is  to  be  referenced  back  to  the 
radiation  input  field  of  the  array.  Thus  the  MTF  is  to  be  modeled  in  terms  of 
the  spatial  frequency  content  of  the  output  of  the  array  compared  to  the  spatial 
frequency  content  at  the  input  surface  of  the  array  itself.  The  NEP  is  to  be 
a  measure  of  that  radiant  input  to  the  focal  plane  surface  which  produces  a 
signal  at  the  output  equal  to  the  noise  level  in  the  output. 

The  thermal  model  described  in  the  preceding  section  models  all  thermal 
processes  and  generates  a  volume  average  temperature  for  use  in  computing  the 
detector  output.  In  particular,  the  thermal  model  computes  a  temperature 
difference  output  signal,  synchronized  with  the  chopper  motion.  This  quantity 
is  the  difference  between  the  volumetric  average  temperature  at  the  time  when 
the  chopper  just  uncovers  the  element  and  at  the  time  it  just  recovers  the 
element.  It  is  computed  for  each  element  throughout  the  thermal  modeling 
process.  When  converted  to  a  voltage,  this  signal  provides  a  rectified  signal  that 
is  synchronized  with  the  chopper  motion.  In  actual  detector  arrays,  this 
rectification  process  is  performed  in  the  array  itself.  This  allows  the  simulation 
of  the  thermal  details  of  the  chopping  process  without  requiring  a  detailed  model 
for  the  rectification  process  of  the  readout. 

These  temperature  differences  can  be  used  along  with  an  overall  model  for 
the  array  to  model  the  performance  of  the  array  in  terms  of  the  electrical 
signals  at  its  output  compared  to  its  radiant  input.  The  concepts  of  MTF, 
responsivity,  and  noise  equivalent  input  all  assume  the  existence  of  a  linear 
relationship  between  the  input  irradiance  and  the  output.  Thus  we  have  enforced 
a  linear  model  throughout  the  analysis. 


Electrical  Modeling 

The  electrical  modeling  consists  of  the  development  of  a  model  whereby 
the  temperature  differences  computed  in  the  thermal  model  can  be  converted  to 
electrical  outputs,  and  the  modeling  necessary  to  properly  introduce  the  noise 
mechanisms  in  the  array  element  and  readout.  The  resulting  model  can  then 
include  the  effects  of  various  noise  sources  on  the  sensitivity  of  the  device  and 
on  its  spatial  resolution. 

The  following  discussion  is  based  on  the  CCD  readout  structure  used  by 
HEOD.  For  modeling  purposes,  the  radiation  pattern  is  stationary  in  space 
except  for  an  on/off  modulation  imposed  by  a  chopper.  The  following 
characteristics  are  carried  over  from  the  HEOD  structure  [9]. 

1)  The  chopper  has  a  50  percent  duty  cycle;  this  is  not  a  restriction  of  the 
thermal  model. 

2)  If  the  chopper  period  is  T  and  the  readout  device  frame  time  is  x,  then 


2Kx  -  T 


where  K  is  a  selectable  integer  value  from  1  to  16. 

3)  The  CCD  readout  device  used  indirect  charge  injection  at  its  parallel 
inputs,  with  an  input  structure  something  like  that  shown  below  in  Figure 
2. 


4)  All  of  the  detector  signals  over  the  whole  focal  plane  are  multiplexed 
through  a  single  output  port  as  described  below.  Figure  3  shows  the 
pattern  of  signal  flow. 
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FIGURE  2.  REPRESENTATIVE  INDIRECT  CHARGE  INJECTION 
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FIGURE  3.  CCD  DATA  FLOW. 


The  readout  cycle  is  as  follows.  The  detector  signals  are  parallel-transferred 
into  the  CCD  structure  of  Figure  3,  so  that  the  MxJ  cells  comprise  an  image 
of  the  detector  signals  integrated  over  a  frame  time.  All  of  the  signals  are 
transferred  to  the  right  by  one  step  (one  cell  width).  Then  the  vertical  CCD 
register  at  the  extreme  right  is  strobed  until  it  has  emptied  all  of  its  data  into 
the  output  port.  Another  single  step  of  the  data  to  the  right  follows,  the 
vertical  CCD  register  is  completely  emptied,  etc.,  until  all  of  the  data  from  the 
focal  plane  has  been  read  out  to  the  CCD  output. 

FOCAL  PLANE  OUTPUT  VOLTAGES 

The  thermal  analysis  code  calculates  a  composite  change  in  detector 
temperature,  AT,  for  each  detector  after  each  time  step, At,  specified  for  the 
finite  difference  thermal  model.  Without  loss  of  generality,  we  can  define  the 
time  t  to  be  zero  at  the  opening  of  the  chopper  (that  is,  the  focal  plane  is  first 
illuminated  by  the  imposed  radiation  pattern  at  t=to=0).  Also  without  loss  of 
generality,  we  can  specify  the  voltage  across  each  of  the  detectors  to  be  zero 
at  t=0.  For  convenience,  and  also  without  loss  of  generality,  the  size  of  a  time 
step  At  can  be  defined  as  an  integer  fraction  of  the  frame  time,  x : 


At  “  —  ,  where  N  is  an  integer.  (18) 


Referring  to  the  circuit  of  Figure  2,  this  same  indirect  charge  injection 
circuit  is  modeled  in  a  form  which  is  appropriate  for  the  necessary  circuit 
analysis  in  Figure  4.  In  Figure  4,  Cj  is  the  sum  of  the  detector  capacitance, 
the  MOSFET  common  source  input  capacitance  and  stray  capacitance.  Rj  is  the 
net  leakage  resistance  shunting  the  detector. 

During  the  n1*1  time  step,  the  (j.m)**1  detector  is  shown  by  the  thermal 
analysis  code  to  have  experienced  a  composite  temperature  change  ATjim  n.  In 
response  to  that  temperature  change,  the  detector  generates  suffficient 


displacement  current  through  its  dielectric  that,  over  the  corresponding  time 
step  of  size  At,  a  charge  of  Aq  has  been  transported  across  the  detector.  Under 
the  reasonable  assumption  that  detector  current  is  a  piecewise  linear  function  of 
time  during  any  time  step,  the  detector  current  (>0)3, m,n  is  defines  as: 

‘  Vl  ‘  i  119) 
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(20) 


where  P  is  the  pyroelectric  coefficient  and  Aq  is  the  area  of  the  detector.  The 
voltage  across  the  (j,m)th  detector  after  the  nth  time  step  is  given  by 


a.  At  - 
j.m.n 


CT  “j.m.n 


(21) 


This  voltage  also  appears  at  the  gate  of  a  common  source  MOSFET  associated 
with  the  (j,m)th  detector.  The  drain  current  of  this  MOSFET  is  related  to  its 
gate  voltage  Fj>m(tn)  by 


•Drain  =  Sm  vj,m,^n^  +  A 


(22) 


where  gm  is  the  MOSFET  transconductance  and  A  is  the  constant  which  is 
related  to  the  MOSFET  threshold  voltage.  Since  the  MOSFET  transconductance 
is  itself  proportional  to  the  square  root  of  iQrajn,  Equation  22  is  nonlinear. 


The  charge  which  is  actually  injected  into  and  transported  through  the 
CCD  is  the  charge  accumulated  on  Cjnt  due  to  idrain*  At  the  beginning  of  each 
frame,  the  reset  MOSFET  is  briefly  strobed  into  the  "on"  state,  so  that  the 
charge  of  cjnt  is  initialized  to  zero.  At  the  end  of  a  frame,  the  discrete  charge 
packet  resulting  from  the  integration  of  i<jrain  over  the  frame  time  is  injected 
into  the  CCD  readout  structure  where  it  subsequently  will  be  transported  to  the 
CCD  output.  A  sample-and-hold  operation  is  typically  performed  at  the  output, 
thereby  transforming  the  charge  packet  into  a  voltage  analog.  The  operation  is 
then  repeated  over  the  span  of  a  subsequent  frame. 

Such  detail  regarding  the  CCD  input  structure  might  be  considered  to  be 
beyond  the  scope  of  the  TFPA  study.  However,  two  important  pieces  of 
information  can  be  extracted  from  this  operational  detail.  The  first  is  that  a 
signal  representative  of  the  integral  of  the  detector  signal  over  a  frame  time  is 
the  quantity  transported  through  the  CCD.  Thus,  the  form  of  the  detector 
outputs  which  must  be  provided  to  the  CCD  is  known.  It  is  simply  the 
checkerboard  image  of  the  focal  plane,  wherein  each  detector  signal  is  updated 
to  the  current  value  at  the  end  of  each  frame.  The  second  is  that  the  indirect 
charge  input  structure  has  some  undesirable  characteristics.  Most  serious  is  the 
nonlinearity  in  Equation  22.  The  concept  of  an  MTF  is  valid  only  for  linear 
systems.  HEOD's  CCD  structure,  after  which  this  exposition  is  patterned, 
incorporates  an  output  compensation  scheme  which  is  intended  to  remove  the 
signal  nonlinearity.  Because  such  a  compensation  scheme  must  operate  upon  the 
signals  after  they  have  been  integrated  over  time,  it  is  clear  that  the 
compensation  scheme  must  inherently  be  imperfect. 

Again,  considering  the  scope  of  the  TFPA  study,  we  have  little  choice  but 
to  assume  that  the  CCD  throughout  is  linear.  The  user  of  the  TFPA  model 
should,  however,  be  aware  of  its  limitations.  When  the  thermal  calculations  are 
performed  at  each  finite-difference  time  step,  Equation  21  can  be  updated  for 
each  detector  with  negligible  additional  computational  load.  After  the  end  of 
a  frame  time,  composed  of  N  time  steps,  has  been  reached,  the  voltage  across 
each  detector  in  the  focal  plane  has  been  determined. 


One  of  the  most  important  characteristics  of  the  array  is  its  MTF,  which 
is  a  measure  of  the  transfer  of  information  through  the  array  as  a  function  of 
spatial  frequency.  Both  thermal  and  electrical  effects  have  an  influence  on  the 
spatial  transfer  function.  Since  we  are  considering  only  rectilinear  arrays  with 
no  thermal  asymmetry  in  this  development,  the  thermal  portions  of  the  MFT  will 
be  identical  along  the  two  dimensions  of  the  array.  Although  it  is  highly 
probable  that  the  two  dimensional  MTF  (wx,Wy)  (where  wx  and  Wy  are  the  spatial 
frequencies  in  the  two  orthogonal  array  dimensions)  is  not  separable  into 
independent  one-dimensional  MTF's,  traditional  modeling  assumes  such  a  sepa¬ 
rability,  and  for  comparison  with  existing  measurements,  we  have  implemented 
only  a  one-dimensional  MTF  (although  the  thermal  code  could  easily  support  a 
two-dimensional  calculation  if  appropriate  sources  were  introduced). 

Consider  the  MTF  for  the  thermal  differences,  AT,  or  its  equivalent,  the 
MTF  of  the  individual  array  voltages.  Referring  to  Figure  3,  let  us  consider  the 
MTF  in  the  horizontal,  or  i^h,  direction.  A  straightforward  way  to  calculate  the 
MTF  of  system  1  in  the  ith  dimension  is  to  apply  a  stationary,  spatially- 
sinusoidal  radiation  pattern  to  the  focal  plane  having  spatial  frequency  wj  in  the 
ilh  dimension  with  complete  uniformity  in  the  jth  dimension.  Because  of  the 
influence  of  the  CCD  readout  device,  the  MTF  exhibits  a  slight  variation  over 
the  focal  plane  and  so  a  specific  row  of  detectors,  say  the  row,  must  be 
chosen  for  analysis.  The  detector  signal  values  obtained  from  Equation  21  at  the 
end  of  the  first  frame  time,  Vj1>m(tn),  must  be  Fourier  transformed  in  space 
along  the  mth  dimension.  As  has  been  indicated  by  the  integrate-and-transfer 
operation  of  the  CCD,  the  detector  signal  form  in  space  is  a  sequence  of 
rectangular  pulses  with  amplitudes  corresponding  to  the  real  detector  signal 
values  in  the  row  and  pulse  width  equal  to  the  geometrical  detector  center- 
to-center  spacing  along  the  mth  dimension. 

The  spatially  transformed  sequence  of  detector  signals  is  Fj(wj).  Similar¬ 
ly,  the  same  sinusoidal  radiation  pattern  except  with  arbitrary  spatial  frequency 
w^  can  be  run  through  the  thermal  model  and  the  resulting  detector  signals  can 
be  spatially  Fourier  transformed  into  F2<W2).  If  the  input  radiation  sinusoids 
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have  the  same  amplitude,  then  the  sequences  of  F(w)  form  the  MTF  for  the  jth 
row  of  detectors  in  the  m-direction.  In  order  to  follow  the  usual  definition  of 
MTF,  the  spatial  frequencies  w^  are  normalized  to  the  Nyquist  spatial  frequency 
along  the  m-direction,  i.e.,  2n  divided  by  the  center-to-center  detector  spacing. 
Thus, 


|F0(0) 


The  function  Fo(0)  must  be  computed  from  the  temperature  output  from  a 
completely  uniform  radiation  field.  However,  the  approach  of  putting  one 
sinusoid  (of  one  frequency)  through  the  thermal  model  at  a  time  is  unnecessarily 
wasteful  of  computer  resources.  A  much  more  reasonable  approach  is  to  apply 
a  radiation  function  to  the  focal  plane  which  contains  all  frequencies.  Any 
waveform  which  contains  a  discontinuity,  such  as  a  step  function,  would  do. 
However,  the  ideal  radiation  input  waveform  is  a  spatial  delta  function,  since 
the  Fourier  transform  of  a  delta  function  consists  of  all  frequencies  equally 
weighted.  Thus,  the  spatial  Fourier  transform  of  the  output  of  the  Jith  row  of 
detector  signals,  resulting  from  a  source  consisting  of  a  one-dimensional  delta 
function,  Fjj(w),  can  be  used  to  formulate  the  MTF  for  the  row  of  detectors 
as  a  continuous  function  of  normalized  frequency. 
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We  can  now  proceed  with  a  more  general,  less  detailed  model  for  the  whole 
detector-readout  system  from  which  expressions  for  the  system  performance,  in 
terms  of  the  results  of  the  thermal  analysis,  can  be  determined. 


GENERALIZED  MODELING 

Despite  the  detailed  discussion  of  the  previous  sections  concerning  the 
structure  of  readout  mechanisms,  chopper-synchronized  differencing  processes, 
and  non-idealities  in  the  readout  systems,  the  wide  differences  between  various 
array  implementations  make  detailed  treatment  of  the  output  processes  of  any 
individual  design  less  than  desirable.  The  resulting  code  would  be  too  limited. 
To  avoid  the  specialization  of  the  code  to  any  particular  readout  mechanism,  we 
have  designed  a  generalized  model  wherein  noise  contributions  can  be  reasonably 
modeled  and  reasonably  general  expressions  for  responsivity,  NEP  and  MTF,  can 
be  derived. 

To  this  end,  we  have  assumed  that  two  linear  system  functions  exist  for 
the  transformation  of  radiant  input  power  from  a  static  (temporally  invariant) 
irradiance  scene  into  the  output  voltages.  (While  the  thermal  model  can  handle 
general  temporally  varying  scenes,  the  present  formulation  for  MTF  is  restricted 
to  time  invariant  patterns.)  The  first  of  these  transfer  functions  is  taken  to 
encompass  all  processes  of  spatial  integration,  pyroelectric  effects,  and  the 
actions  of  the  chopper,  including  the  signal  detection  process  in  synchronism 
with  the  chopper.  Noise  contributions  include  all  those  noises  up  to  the 
extraction  of  a  sample  voltage  level.  Thus,  the  output  of  this  first  process 
consists  of  an  array  of  voltage  sequences,  a  three-dimensional  function,  two 
dimensions  of  which  are  the  array  spatial  dimensions  and  the  third  of  which  is 
the  sample  sequence  dimension,  corresponding  to  time.  Note  that  in  the 
generation  of  these  sequences,  and  in  the  modeling  of  the  effects  of  the  detector 
processes,  no  reference  to  time  is  needed;  the  detector  outputs  are  simply 
discrete  voltage  outputs.  The  noises  which  contribute  to  this  voltage  sequence 
consist  of  all  noise  processes  present  in  the  array  elements  themselves.  Since 
the  detector  elemental  voltages  are  sampled  to  obtain  the  output  voltage 
sequences,  the  expressions  for  the  power  spectra  of  the  various  noise 
contributions  must  be  modified  to  account  for  the  integration  process,  and  the 
variance  of  the  sampled  noise  computed. 

The  resulting  sequences  of  signal  voltages  (and  noise  voltages)  are  then 
passed  through  the  second  process,  the  readout  device,  which  orders  the  array 
of  voltage  sequences  into  a  single  sequence  of  voltages  for  passage  out  through 


a  single  line.  This  process  both  distorts  the  signals  and  adds  noise.  These 
effects  of  the  readout  device  contribute  to  both  the  degradation  of  spatial 
resolution  and  to  the  noise  of  the  system. 

The  detector  array  signals  form  a  matrix  in  space  and  a  sequence  in  time, 
with  the  time  duration  of  each  element  of  the  sequence  chosen  so  that  the  whole 
array  is  read  out  in  one  frame  period.  The  sequence  of  voltages  is  not  itself 
examined  in  time,  but  is  presumably  used  to  generate  an  image  of  the  array 
output  for  presentation  to  the  user.  Thus  the  spatial  frequency  content  as 
apparent  in  the  display  of  the  array  output  is  to  be  compared  to  the  spatial 
frequency  content  of  the  radiant  input.  The  time  domain  (for  stationary  scenes) 
plays  only  a  very  incidental  role  in  such  a  process.  Time  is  important  only  to 
the  extent  that  it  takes  time  to  sample  the  array,  order  the  array  of  voltages 
into  a  single  line  and  present  it  to  this  display.  Time  is  also  important  in  that 
the  noise  processes  of  the  detector  elements  and  readout  devices  are  specified 
in  terms  of  temporal  power  spectral,  which  must  be  integrated  to  obtain  sample 
statistics  from  temporal  statistics. 

Thus  an  overall  model  for  the  thermal  focal  plane  array  suitable  for 
definition  of  typical  performance  parameters  can  be  constructed  as  the  following 
set  of  linear  operations: 

IRRADIANCE  PATTERN  IN  SPACE 

DETECTOR  NOISE 

TWO-DIMENSIONAL  ARRAY  OF  VOLTAGE  DIFFERENCES 

READOUT  NOISE 

SINGLE  SEQUENCE  OF  VOLTAGE  VALUES 

* 

ARRAY  OF  BRIGHTNESS  VALUES  ON  A  DISPLAY 


The  vertical  arrows  represent  the  linear  processes;  the  inserted  noises  represent 
the  origination  of  noise  contributions.  Two  noise  contributions  are  important: 
those  which  originate  in  the  detector  elements  themselves  and  those  which 
originate  in  the  readout  device.  Johnson  noise,  so-called  "loss  tangent"  noise  and 
noise  due  to  thermal  fluctuations  are  typical  noises  which  originate  in  the 
detector.  (Individual  noises  are  considered  in  a  later  section.)  Amplifier  noise, 
reset  noise  and  other  electronic  uncertainties  constitute  the  contribution  labelled 
readout  noise.  In  the  analysis  presented  here,  all  noise  processes  have  been 
referenced  to  the  readout  input.  Thus,  all  noises  can  be  lumped  together  and 
assumed  to  originate  at  the  detector.  This  means  that  the  processes 
representing  the  readout  and  display  processes  can  also  be  lumped  together. 
Figure  5  shows  a  "black  box"  flowchart  of  these  processes. 

In  this  analysis,  we  have  chosen  to  ignore  all  temporal  variations  in  the 
input  scene  and  assume  the  scene  to  be  static.  Our  neglect  of  these  variations 
is  prompted  by  the  required  chopper  system;  it  greatly  complicates  the 
treatment  of  temporally  varying  scenes. 

The  scene  is  considered  static  for  this  performance  parameters  (MTF  and 
NEP)  evaluation,  so  that  only  spatial  variations  are  of  interest.  Within  this 
limitation  and  in  the  spirit  of  other  definitions,  the  MTF  of  the  system  will  be 
defined  as  the  filtering  performed  by  the  overall  system  on  spatial  information 
expressed  in  two-dimensional  spatial  frequency  space. 

NEP  or  NEFD  or  equivalently,  noise  equivalent  irradiance  is  an  expression 
of  the  radiant  input  required  to  give  an  output  equivalent  to  the  noise  generated 
within  the  device.  It  can  also  be  thought  of  as  a  function  in  two-dimensional 
spatial  frequency  space.  The  noise  must  be  considered  as  the  noise-generated 
frame-to-frame  variations  in  the  display  device  output. 

Finally,  the  responsivity  is  a  measure  of  the  conversion  of  radiant  energy 
into  electrical  signal,  generally  having  units  of  volts/watt  or  volts/watt/cm2. 
This  can  also  be  thought  of  as  a  spatial  frequency  dependent  quantity  and  is 
generally  equivalent  to  the  transfer  function  of  an  individual  detector  element, 
i.e.,  the  transfer  function  of  the  leftmost  box  in  Figure  5.  The  commonly  quoted 
value  is  the  DC  or  zero  spatial  frequency  value. 
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Each  of  the  functions  represented  by  boxes  in  Figure  5  is  assumed  to  be 
a  linear,  spatially  invariant  process  so  that  the  actual  output  of  each  box  can 
be  represented  as  a  convolution  of  its  input  with  its  transfer  function  in  space. 
When  this  is  the  case,  the  output  of  each  box  in  two-dimensional  frequency  space 
can  be  represented  as  the  product  of  the  (two-dimensional)  Fourier  transform  of 
its  input  and  the  (two-dimensional)  Fourier  transform  of  its  transfer  function. 
The  transformation  between  the  physical  space  of  the  input  scene  and  the  two- 
dimensional  array  index  space  of  the  array  is  included  in  the  convolution  process. 
Similarly  the  transformation  from  the  two-dimensional  array  space  of  the 
detector  element  to  the  one-dimensional  output  sequence  and  back  to  spatial 
display  space  is  represented  by  a  convolution  process.  The  corresponding  Fourier 
transforms  give  the  relationships  between  the  various  frequency  spaces. 

The  following  equations  develop  these  relationships  in  more  detail  in  order 
that  precise  expressions  for  MTF,  NEFD  and  responsivity  can  be  developed  in 
terms  of  the  thermal  code  outputs,  the  array  physical  properties  and  the 
postulated  readout  characteristics. 

The  two-dimensional  array  of  detector  output  voltages  can  be  determined 
from  the  irradianee  input  image  as  follows: 


V(i,j) 


LL 


Ha(1Ax-x' , jAy-y')  E(x',y')  dx'dy’ 
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Here  is  linear  process  transducing  the  irradianee  distribution  E(x,y)  incident 
on  the  array  into  the  i,j  array  of  detector  voltages.  Note  that  includes  all 
chopper  effects  and  all  effects  of  detection  area,  x  and  y  are  the  dimensions 
of  an  individual  array  element.  Appendix  C  presents  an  initial  description  of  the 
system  including  the  details  of  the  chopper  operation. 

This  voltage  signal  is  passed  through  a  readout  device  and  is  then  displayed 
as  a  two-dimensional  display  in  which  each  element  of  the  two-dimensional  array 
of  voltages  is  converted  to  an  output  value  corresponding  to  the  appropriate 
pixel.  Thus  although  the  image  output  consists  of  a  discrete  array  of  voltages, 
the  frequency  information  can  be  kept  in  the  spatial  domain. 


Furthermore,  a  continuous  representation  of  detector  voltages  is  preferable 
to  the  discrete  representation.  It  avoids  the  problems  associated  with  the 
conversion  to  and  from  the  discrete  space  of  the  detector  elements.  The  voltage 
as  a  continuous  function  of  spatial  position  (Equation  25)  can  then  be  written  in 
the  spatial  domain  as 


V(x,y) 


If 


HA(x-x' ,y-y’ )  E(x',y')  dx'dy' 


(26) 


In  the  spatial  frequency  domain,  this  becomes 

v<wx»wy>  "  WV  *  E(VV  (27) 

where  functional  dependence  on  a  pair  of  frequencies  implies  the  application  of 
a  Fourier  transform  to  the  function.  wx  and  Wy  are  spatial  frequencies  in  the 
x-  and  y-  directions  respectively.  Thus 
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and  similarly 
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H^(wxWy)  is  the  spatial  frequency  dependent  responsivity  function  for  the  focal 
plane  array  alone,  i.e.,  exclusive  of  the  effects  of  the  readout  device.  The  so- 
called  DC  responsivity  of  the  array  is  H ^(0,0)  and  the  MTF  of  the  array  is  the 
amplitude  of  the  response  at  any  given  frequency  compared  to  that  at  DC. 
Thus,  the  magnitude  of  the  responsivity,  is: 

RA(VV  “lHA(VV'  (31) 
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The  only  effects  of  the  readout  device  on  the  MTF  that  we  have  included 
in  the  present  model  are  those  of  incomplete  charge  transfer  between  CCD 
cells.  These  effects  are  properly  considered  in  a  one-dimensional  discrete 
domain.  However,  for  the  CCD  systems  to  be  used  with  these  arrays,  the 
transfer  arrangement  is  such  that  all  transfers  in  one  direction  are  made  before 
any  in  the  other.  Thus  the  smearing  effects  of  charge  transfer  remain 
completely  separated  in  the  two  dimensions,  and  the  effects  of  a  CCD  on 
responsivity  and  MTF  can  be  treated  as  simply  the  product  of  two  independent 
effects,  so  that  the  overall  MTF  can  be  written: 


MTF(v  ,w  )  »  MTF  (w  ,w  )  •  MTF___(w  )  MTF(w  ) 
x  y  A  x  y  CCD  x  CCD  y 


The  two  CCD  effects  can  be  written  [10] 


1,11  CCD (Wx>  “  exp  ^‘xli  “  eo»(wAx)]} 


MTFr™<'0  *  «*P  11  -  cos(v..  y)]} 


where  1  and  J  are  the  number  of  transfers  in  the  x-  and  y-directions  respectively 
and  ex  and  e  y  are  the  respective  single  stage  transfer  losses.  Note  that  since 
MTFCCD«»  =  1  for  values  of  le*,  the  charge  transfer  inefficiencies  do  not 
contribute  to  the  DC  responsivity.  Thus 

R(0,0)  -  R. (0,0)  <37) 
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H^Wx.Wy)  is  obtained  as  the  Fourier  transform  of  H^USy’),  which  is  in 
turn  determined  from  the  output  of  the  thermal  code.  First  of  all,  to  make  the 
results  of  the  code  comparable  to  other  modeling  studies,  we  have  made  the 
common  assumption  that  H^(wx,Wy)  is  separable  into  two  independent  one¬ 
dimensional  contributions  for  the  x  and  y  coordinates  respectively.  By  making 
this  assumption,  we  can  use  a  line  source  to  represent  a  delta  function  source 
in  one  dimension  and  use  it  to  probe  the  response  function  of  an  array  element 
in  that  direction.  Note  that  an  integral  over  a  delta  function  can  be  viewed  as 
"pulling  out"  the  value  of  the  remainder  of  the  integrand  at  the  point  where  the 
delta  function  is  applied.  Thus  the  use  of  a  delta  function  in  irradiance  as  input 
to  the  thermal  code  can  be  used  to  determine  values  of  H^(x)  as  a  function  of 
x.  Under  the  separability  assumption 

wv  ■  ha  <".>  aA  <v  <39) 

x  y  J 

~nd  with  the  use  of  a  line  source  in  the  y-direction  at  the  position  xQ, 
i.e., 

E(x,y)  -  Ey  5(x-xq)  (40) 
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where  Ey  is  the  irradiance  of  the  line  source  per  unit  source  length,  we  obtain 
from  Equation  25, 
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or  in  terms  of  a  separable  H: 
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(Index  j  becomes  irrelevant  here  since  all  elements  j  will  see  identical  results.) 
The  output  of  the  thermal  code  is  used  to  determine  V(i)  for  various  positions, 
x0,  of  the  input  delta  function  source.  Then  can  be  determined  from 
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or  by  simply  changing  variables, 
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Note  that  the  integral  of  the  response  in  the  y-direction  over  the  whole  y- 
direction  is  a  constant  in  this  determination  of  the  x-direction  response. 

A  determination  of  the  response  in  the  y-direetion  can  be  performed 
similarly  but  with  the  line  source  in  the  x-direction.  With  present  thermal  code, 
this  is  not  generally  done  since  it  is  simpler  to  turn  the  array  and  perform 
another  set  of  calculations.  Furthermore,  if  the  array  elements  are  square, 
there  is  no  need  to  run  a  second  direction;  the  isotropy  of  the  materials  ensures 
identical  results. 

The  DC  responsivity  of  the  system  is  simply  the  integral  of  H^(x,y)  over 
all  space  x  and  y.  Within  the  separability  assumption,  the  responsivity  at  DC 
can  be  determined  from  Equation  45  as. 


(46) 


The  responsivity  is  simply  the  integral  of  the  output  voltage  translated  into 
space  divided  by  the  line  source  intensity  in  watts/cm.  Note  that  the 
responsivity  can  be  obtained  from  results  in  only  one  direction. 

Furthermore,  the  MTF  function  can  be  obtained  from  a  Fourier  transform 
of  V(x),  normalized  to  its  zero  frequency  value. 


MTF  (v  ) 

X  X 
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To  obtain  an  MTF  in  one  direction,  it  is  not  necessary  to  obtain  a  value  for  the 
integrated  response  function  in  the  other  direction. 


Formulation  of  Noise  Components,  NEP  and  D*  for  Pyroelectric 

Detector/CCD  Readout 

Such  important  system  parameters  as  NEP  and  D*  depend  upon  appropriate 
formulations  of  system  noise.  The  following  is  the  formulation  of  relevant  noise 
sources  as  well  as  formulation  of  NEP  and  D*  expressions.  The  noise  equivalent 
model  of  the  system  is  shown  in  Figure  6. 

Noise  is  introduced  into  this  system  by  the  sampling  of  the  time  dependent 
noise  generated  both  in  the  detector  elements  and  in  the  input  stages  of  the 
system  electronics.  The  sampling  of  this  noise  leads  to  uncertainties  in  the 
output  voltage  array  which  vary  from  frame- to- frame.  Thus,  the  noise  in  the 
system  can  be  characterized  by  the  temporal  noise  power  spectrum  of  each  of 
the  noise  contributions,  which  can  be  added  to  yield  a  total  noise  spectrum.  The 
composite  noise  as  it  appears  in  the  system  output  is  characterized  by  this 
composite  noise  spectrum  as  integrated  over  a  sample  time.  The  noise  passes 
only  through  the  readout  device;  therefore  the  only  frequency  content  modifica¬ 
tions  are  made  by  that  device. 

We  assume  that  all  noises  can  be  referenced  to  the  detector  output.  HEOD 
[9]  has  advised  us  that  the  MOSFET  used  for  indirect  charge  injection  into  their 
CCD  has  sufficient  gain  that  noise  sources  downstream  of  the  MOSFET  can  be 
ignored.  Similarly,  we  assume  that  noise  sources  downstream  of  the  T1 
multiplexer  can  wither  be  referenced  to  the  multiplexer  input  or  ignored.  Thus, 
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FIGURE  6.  SYSTEM  NOISE  EQUIVALENT  MODEL. 
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the  noise  current  power  spectra  can  be  summed  in  quadrature  into  a  single  noise 
power  spectrum.  So 


1(f)  - 


l/(f) 


(48) 


and  the  components  of  1(f)  are  defined  to  be: 
=  Johnson  noise  current  of  R 

(  4kT  V*  amp 

\  *4/  /Sz 

12  =  radiation  noise  current 
-  (16  kor^)4^  ^ 


I3  =  Johnson  noise  of  detector  dielectric 


■  —  (BnkTA.e'e  f  tan  6)^ 

/d  d  °  /Hz' 
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I4  =  Statistical  noise  due  to  discrete  nature  of  gas  conduction  between  detectors 

■  MW*  “ 


I5  =  Statistical  noise  due  to  lateral  thermal  conduction 
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Ig  =  Equivalent  input  noise  current  of  indirect  charge  injection  input 
MOSFET 


specify  default^  5xl0-17 

/Hz 


I7  =  The  equivalent  of  input  noise  voltage  of  input  MOSFET 


where:  k  =  Boltzmann's  constant  =  1.38xl0"22  (watts-see/K) 

0  =  Stefan-Boltzmann  constant  =  5.67x10“®  (watts/m^-K^) 

T  =  Detector  nominal  temperature  (Kelvin) 

Ry  =  Net  leakage  resistance  across  detector 
R-p  =  User  specify  default^  jq13  (^) 

A(j  =  Detector  area  (m2) 
a  =  Absorptivity  of  detector  active  surface 
P  =  Pyroelectric  constant  (Q/m2-K) 

C’  =  volumetric  specific  heat  (J/m®-K) 
d  =  Detector  thickness  (m) 

e'  =  Real  part  of  detector  relative  dielectric  constant 

=  Free  space  dielectric  constant  =  0.111x10"®  (Farad/m) 

f  =  electric  frequency  (Hz) 

tan  6  =  Dielectric  loss  factor 

K  =  Thermal  conductivity  of  gas  surrounding  detectors 
(watt/m-K) 

kg  =  Thermal  diffusivity  of  gas  surrounding  detectors  (m2/s) 
Kp  =  Thermal  conductivity  of  detector  material  (watt/m-K) 


kp  =  Thermal  diffusivity  of  detector  material  (m2/s) 

|Z|  =  Steady-state  value  of  parallel  impedence  across  detector 

(G) 
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where:  C-p  =  Total  parallel  capacitance  across  detector  (Farad) 

=  CD  +  Cs 

Cp  =  Detector  capacitance  (Farad) 

e’e  A. 

o  d 


Cs  =  Total  parallel  stray  capacitance  (Farad) 
=  User  specify  default^  0.15xl(Tl2  p^ad 


B  *  10 


(  30  «,  v 

v° + V 


where:  f0  =  Corner  frequency  of  input  MOSFET  1/f1/2  noise 
=  User  specify  default^  an  Hz 


The  noise  generated  at  the  input  of  the  CCD  is  that  uncertainty  resulting 
from  the  integration  and  sampling  of  the  noise  spectrum,  resulting  in  a  single 
random  quantity  for  each  voltage  measurement.  We  can  reasonably  assume  that 
these  sampled  noise  voltages  will  be  uncorrelated  if  the  array  has  a  reasonably 
large  number  of  elements.  When  the  number  of  elements  is  large,  the  time 
between  each  sample  from  any  one  element  must  then  be  larger  than  a  single 
sample  time,  and  if  the  appropriate  filtering  has  been  performed  in  the  input 
stages  to  eliminate  having  time  correlations  much  larger  than  a  sample  time, 
then  the  samples  will  indeed  be  uncorrelated.  Of  course,  this  assumption 
requires  the  presence  of  such  low  pass  filtering  and  subsequent  sections  include 
such  filtering,  but  we  will  ignore  correlations  between  noise  samples  in  the 
analysis  that  follows.  These  uncorrelated  noise  samples  are  associated  with 


i 

i 


35 


individual  elements  of  the  output  voltage  array,  and  the  spatial  power  spectrum 
of  the  noise  input  to  the  readout  structure  is  simply  the  variance  of  these  noise 
currents  samples,  multiplied  by  the  square  of  the  parallel  input  impedence  across 
the  detector,  and  by  the  sine  squared  x  over  x  squared  spatial  dependence 
determined  by  the  rectangular  array  elements.  We  obtain  the  following  form  for 
the  noise  spatial  power  spectrum  as  seen  in  the  output  device: 
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where  oj^(At)  is  the  variance  of  a  composite  noise  current  sample  taken  over  time 
At.  The  determination  of  0j2(At)  will  be  discussed  subsequently.  The  factor  of 
2  arises  because  each  output  value  results  from  two  voltage  measurements,  one 
from  each  state  of  the  chopper. 

The  effect  of  the  readout  structure  is  to  modify  the  noise  power  spectrum 
in  spatial  frequency  by  the  MTF  of  the  CCD  structure.  Thus, 
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The  NEP  of  the  system  is  defined  as  that  radiant  input  required  to  produce 
an  output  amplitude  just  equal  to  the  noise  amplitude  at  the  system  output. 
Equation  60  gives  the  output  in  terms  of  the  radiant  input.  Equating  the  output 
from  a  source  NEP(wx,Wy)  to  the  noise  output  yields  the  appropriate  expression 
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This  is  the  amplitude  spectrum  of  the  noise  equivalent  power  in  the  spatial 
frequency  domain,  and  represents  the  radiant  input  at  any  spatial  frequency 
required  to  reproduce  the  noise  amplitude  apparent  between  individual  frames, 
at  that  spatial  frequency.  Thus  the  temporal  bandwidth  of  this  NEP  spectrum 
will  be  that  associated  with  the  frame  time  of  the  array.  The  total  NEP  is  just 


the  DC  radiant  input  required  to  produce  the  (spatial)  variance  in  output,  and  is 
thus  the  integral  of  the  noise  spectrum  over  all  spatial  frequencies  divided  by 
the  DC  responsivity,  thus 


1/2 


1/2 


(63) 

Note  that  the  spatially  dependent  NEP  is  not  constant  with  spatial 
frequency.  Thus,  the  NEP  will  not  generally  be  the  same  if  the  power  were  to 
be  applied,  say,  over  a  square  pattern  and  over  a  long,  narrow  rectangular 
pattern.  It  would  not  even  be  the  same  if  the  power  were  to  be  applied  over 
two  square  patterns  of  different  sizes.  Such  quantities  as  noise-equivalent 
temperature  (NEtfT)  and  minimum  resolvable  temperature  (MRT)  would  similarly 
be  associated  with  specific  spatial  patterns  of  radiation. 

Values  for  D*  can  be  computed  for  the  system  as  follows.  To  the  extent 
that  the  thermal  focal  plane  array  has  absorptivity  which  is  spectrally  uniform, 
then  D*biac|<body  =  Dp*  =  D*  is  formulated  as 
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Finally,  it  is  necessary  to  determine  0j2(^t).  The  noise  current  spectrum 
is  the  sum  of  the  noise  current  spectra  of  each  of  the  contributions,  Ij,  described 
earlier.  Since  the  noise  contributions  all  have  a  mean  of  zero,  the  mean  value 
of  the  noise  sample  will  also  be  zero.  But  the  variance  of  the  noise  sample  will 
depend  on  the  temporal  noise  spectrum  and  the  time  over  which  the  integral  is 
taken.  A  value  for  the  variance  can  be  determined  by  noting  that  integration 
over  a  time  interval  At  is  equivalent  to  integration  over  a  unit  step  function  of 
duration  At 
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where  Mt)  is  a  unit  step  of  duration  At,  i.e.. 
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But  the  effect  of  this  step  function  on  the  power  spectrum  is  to  multiply  it  by 
[(sin  wAt/2)/(wAt/2)l2.  Thus 


Sy  (W) 
int 


S„  (w) 


/sin 

\  0 


(wAt/w)  \ ' 
(wAt/2)  J 


(66) 


Furthermore,  the  variance  of  a  function  is  just  the  integral  of  its  power 
spectrum  over  all  frequencies.  In  terms  of  the  summed  noise  spectra  1(f) 


o^Wt) 


IW(w/2ir)  /sin  wAt/2 
l  wAt/2 


dw 


(67) 


(The  I,(f)  functions  are  all  single-sided  power  spectra.)  To  determine  Oi2(At)  in 
the  code,  first  a  set  of  w  values  are  selected  based  on  At.  They  have  been 
chosen  to  properly  cover  the  spectral  range  between  DC  and  4ir/At. 

In  practice,  we  have  incorporated  an  additional  low  pass  filtering  operation 
to  eliminate  all  noise  substantially  below  the  chopping  frequency  of  the  system. 
This  is  accomplished  by  the  multiplication  of  1(f)  by  a  filter  having  a  nominal 
cutoff  at  one-half  the  chopping  frequency.  Thus 


I2(w/2ti) 


(68) 


where  Wc  *s  the  chopping  frequency.  This  low  pass  filtering  operation  will 
prevent  very  low  frequency  noise  from  contributing  to  Oj2  when  in  practice  it 
would  be  filtered  out  of  the  system.  It  will  also  help  to  insure  that  adjacent 
samples  from  any  one  detector  will  be  uncorrelated. 

Note  that  the  use  of  the  quantity  0j2(^t)  as  representing  the  detector  noise 
contribution  signifies  the  transformation  of  the  electrical  noise  generated  within 
the  detector  element  (and  of  other  input  processes  referenced  back  to  the  input) 
into  discrete  voltage  samples  which  are  then  passed  through  the  discrete  CCD 
svstem  and  a  discrete  element  output.  Thus  the  discrete  independent  noise 
values  are  smeared  by  the  CCD  operations  and  mapped  onto  the  display.  The 
spatial  spectrum  of  the  independent  noise  samples  is  now  given  by  the  spatial 
spectrum  of  the  array  (or  equivalently  of  the  display)  as  modified  by  passage 
through  the  CCD. 

In  summary,  the  electrical  portion  of  the  code  acts  like  a  postprocessor, 
using  a  series  of  line  sources  at  several  positions  to  probe  the  lumped  impulse 
response  function  of  an  individual  detector  element.  The  MTF  of  the  array  is 
determined  by  the  spatial  Fourier  transform  of  this  sampled  impulse  response. 
The  responsivitv  is  determined  from  the  integration  of  this  impulse  response  over 
the  probed  dimension.  The  noise  of  the  system  is  modeled  using  a  variety  of 
contributions  to  detector  element  noises,  plus  the  input  noises  referenced  to  the 
detector  element  outputs. 

Finally,  the  spatial  frequency  dependent  NEP  and  D*  values  are  computed 
from  the  responsivity,  MTF  and  a  sampled  noise  variance  value  computed  using 
the  composite  detector  noise  spectrum.  The  total  NEP  is  based  on  the  signal 
necessary  to  produce  the  total  output  variance. 
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SECTION  5.  CODE  DESCRIPTION 


The  thermal  focal  plane  array  code  has  three  main  parts: 

1.  The  Preprocessor  produces  an  efficient  set  of  input  for  the  thermal 
calculations  given  the  array  configurations,  dimensions  and  materials. 

2.  The  Processor  calculates  the  transient  three-dimensional  temperature 
distribution  in  array  given  inputs  from  preprocessor  and  boundary 
condition  specification  (radiative,  convective  and  conductive).  It  also 
calculates  voltage  outputs  of  all  detectors  as  a  function  of  time  given 
their  electrical  and  pyroelectric  properties. 

3.  The  Postprocessor  evaluates  various  detector  performance  parameters 
from  the  calculated  voltage  responses  (generally  for  specific  boundary 
conditions  and  configurations  treated  in  the  processor  and  preprocessor 
respectively)  considering  various  thermal,  electrical  and  signal  pro¬ 
cessing  noise  sources. 

Each  of  these  components  is  discussed  in  order  in  the  following  three 
subsections.  The  numerical  procedures  and  approximations  will  be  discussed  in 
detail. 


5.1  PREPROCESSOR 

The  only  purpose  of  the  thermal  focal  plane  array  preprocessor  is  to  set 
up  input  files  for  the  processor  and  postprocessor.  This  module  removes  most 
of  the  effort  in  setting  up  the  finite  difference  node  structure  and  heat  transfer 
paths,  and  optimizes  the  files  for  maximum  computational  efficiency  in  the 
processor. 

The  governing  differential  equation  (Equation  1)  is  expanded  in  finite 
difference  form  [6]  to  give  for  an  arbitrary  node: 

n 

PC  V[T(t+At)  -  T(t) ]  -  At  l  Q  (69) 

P  i-1  1 

where  pCpV  is  the  volumetric  heat  capacity  of  the  node  (e.g.,  J/K);  T  is  the  node 
temperature  at  discrete  time  points  (e.g.,  t,  t+At,  t+2$,  ...);  and  Qj  are  the  net 
heat  transfer  rate  to  the  node  (e.g.,  J/s)  from  various  sources  (i  =  1,2,3,...). 
These  heat  transfer  terms  may  be  due  to  only  heat  conduction  (e.g.,  for  internal 


nodes),  or  can  be  due  to  combinations  of  conduction,  radiation  and  convection 
(e.g.,  for  surface  nodes).  The  conduction  terms  are  expressed  simply  as  finite 
difference  in  space: 
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where  k  is  the  thermal  conductivity  of  the  material  between  the  i  and  j  nodes, 
A  is  the  cross-sectional  area  of  the  conduction  path  and  Axjj  is  its  length. 

Note  that  this  expression  evaluates  the  heat  transfer  based  on  the 
temperature  gradient  at  the  start  of  the  time  interval  (i.e.,  at  t).  This 
approximation  introduces  the  usual  finite  difference  errors  and  instabilities,  but 
for  the  check  cases  we  considered  in  the  validation  (see  Section  6),  these  do  not 
appear  to  cause  problems  for  the  temporal  and  spatial  computational  step  sizes 
used.  Other  implicit  techniques  do  not  have  these  potential  problems,  but 
require  more  complex  numerical  procedures  and  so  were  not  used  in  this  initial 
version  of  the  simulation  code. 

The  preprocessor  automatically  subdivides  the  detector  elements  and 
surrounding  materials  into  a  three-dimensional  orthogonal  grid  of  nodes  whose 
temperatures  are  to  be  calculated  by  the  processor.  These  nodes  are  regularly 
spaced  in  each  material,  but  arbitrary  fineness  of  nodes  can  be  specified  for 
each  material  and  each  direction  within  the  practical  storage  and  run  time 
limitations  of  the  computer.  Note  that  spatial  resolution  in  the  calculation 
influences  temporal  step  size  because  of  calculational  stability  considerations. 

The  preprocessor  also  automatically  evaluates  the  effective  path  con¬ 
ductance  for  all  paths  connecting  all  nodes: 


Kjj  —  kA/ Akjj 


(71) 


The  effective  value  calculated  for  paths  between  different  materials 
considers  the  different  material  properties  but  neglects  any  contact  resistance 
between  the  materials.  The  code  than  produces  files  for  the  processor 
containing  the  effective  conductances  for  each  path,  the  identification  of  the 
node  pairs  associated  with  each  path,  and  the  total  heat  capacity  of  each  node: 

Cj  =  pCpV  (72) 

Actually,  the  files  contain  "pointers"  to  efficiently  specify  the  conductance  K 
and  capacity  C.  This  is  done  because  many  of  the  paths  will  have  the  same 
value  of  K,  and  many  of  the  nodes  will  have  the  same  value  of  C.  For  example, 
all  interior  nodes  in  the  detector  may  have  a  capacitance  value  C(n)  while  all 
interior  nodes  in  the  insulation  in  some  region  may  have  a  value  C(m).  Then  the 
files  only  need  to  specify  the  indices  or  pointers,  n,  m,  etc.,  along  with  tables 
of  C(j),  etc. 

The  radiative  heat  rates  to  all  the  irradiated  surface  nodes  is  ap¬ 
proximately  by: 

«r.d  ‘  ♦  «<=(*„*  -  V)  1731 

where  a  is  the  absorptance  of  the  surface,  A  is  the  surface  area  of  the  node, 
E  is  the  irradiance  on  the  node  (e.g.,  J*m-2-s~l)  due  to  the  focused  infrared 
radiation  that  is  to  be  measured,  e  is  the  emittance  of  the  surface,  a  is  the 
Steffan-Boltzmann  constant,  and  T5  is  the  effective  temperature  of  the 
radiative  environment  including  everything  except  the  focused  radiation.  This 
assumes  that  all  the  surface  nodes  exchange  radiative  energy  with  the  same 
constant  and  uniform  background. 

The  preprocessor  automatically  specifies  the  nodes  that  are  on  this 
irradiated  surface,  and  calculates  tables  of  radiative  exchange  factors  for  each 
node  type  (e.g.,  detector  node,  insulator  nodes): 

Fr  =  ctA  Fe  =  eAo  (74) 

These  values  are  again  specified  by  pointers  for  the  sake  of  storage  efficiency. 


The  other  exterior  surfaces  of  the  detector  element  that  are  not  irradiated 
by  focused  radiation  may  also  exchange  heat  with  their  environment  by 
radiation.  For  example,  some  configurations  may  have  the  detectors  fully 
reticulated,  but  not  backfilled  with  insulation,  so  nodes  on  adjacent  detectors 
may  transfer  heat  by  radiation.  Exchange  factors  for  these  nodes  are 
automatically  calculated  by  the  preprocessor  based  on  Equation  11: 

Fs  =  AiF12  (75) 

This  treatment  assumes  that  the  net  radiative  heat  transfer  to  side  surface 
nodes  (i.e.,  not  irradiated  by  focused  radiation)  is  influenced  only  by  the  directly 
opposing  node  on  the  adjacent  detector.  The  net  radiation  from  all  other  sources 
is  assumed  to  be  zero;  that  is,  on  the  average  all  other  surfaces  supply  as  much 
energy  as  they  absorb  (so-called  "refractory"  surfaces).  The  preprocessor 
automatically  specifies  the  radiatively  interacting  node  pairs,  either  between 
directly  opposing  nodes  on  adjacent  detectors,  or  between  nodes  on  the  "bottom" 
surface  of  the  detector  and  the  substrate. 

Convective  heat  transfer  to  surface  nodes  is  also  handled  in  a  similar 
fashion.  Using  Equation  8,  the  code  calculates  convective  exchange  factors  for 
all  pertinent  nodes. 

Fc  =  Ah  (76) 

This  treatment  assumes  the  ambient  gas  is  an  infinite  heat  source  since  its 
temperature  never  changes.  The  correlation  used  as  input  is  based  on 
measurements  in  air,  but  any  convective  coefficient  can  be  specified  as  input. 

The  preprocessor  also  calculates  a  table  of  surface  node  positions. 
Currently,  each  position  is  given  in  terms  of  one  coordinate  only  since  this  is 
how  the  irradiance  function  E  is  currently  specified  in  the  processor.  Extension 
to  two-dimensional  irradiance  patterns  would  require  two  coordinates  to  be 
transferred  for  each  irradiated  surface  node. 


Other  output  files  are  also  prepared  for  the  postprocessor.  These  files 
contain  the  configuration  of  the  array,  the  area  of  the  detectors  and  the  position 
of  the  irradiated  detector. 

5.2  PROCESSOR 

The  processor  reads  the  inputs  from  the  preprocessor  along  with  the 
parameters  that  specify  the  boundary  conditions,  and  calculates  the  temperature 
of  all  the  nodes  at  each  time  step.  This  calculation  is  carried  out  by  simply 
evaluating  the  heat  transfer  for  each  path.  The  resulting  temperature  changes 
of  the  nodes  for  each  path  are  accumulated,  and  after  all  the  paths  have  been 
considered  the  new  temperature  of  each  node  is  evaluated  from  their  initial 
temperature  and  net  temperature  change.  These  calculations  are  performed  for 
as  many  time  steps  as  desired  within  the  computer  run  time  limitations. 

The  processor  also  evaluates  the  volumetric  average  temperature  of  all  the 
nodes  in  each  detector  at  specified  time  increments.  Since  all  the  nodes  in  the 
detectors  have  the  same  volume,  the  volumetric  average  is  simply  the  arithmetic 
average  of  all  the  node  temperatures.  The  voltage  output  of  the  detectors  is 
directly  related  to  the  difference  in  the  volumetric  average  temperature  at  the 
beginning  and  end  of  the  chopping  cycle.  These  voltages  as  a  function  of  time 
for  each  detector  are  stored  in  a  file  for  further  analysis  in  the  postprocessor. 

5.3  POSTPROCESSOR 

The  purpose  of  the  postprocessor  is  to  simulate  the  readout  and  signal 
processing  functions  of  the  thermal  focal  plane  array  system,  and  to  evaluate 
the  performance  of  the  system  in  terms  of  various  criteria.  All  the  equations 
for  the  noise  terms  (see  Section  4.2)  are  straightforward  algebraic  expressions 
that  present  no  problem  for  computer  evaluation.  Similarly  the  expressions  for 
the  performance  parameters  are  also  straightforward.  The  only  problem  comes 
in  evaluating  the  parameters  at  low  chopping  frequencies  where  the  thermal 
solution  requires  large  amounts  of  computer  time.  However,  no  special 
numerical  techniques  are  required  in  the  postprocessor.  The  postprocessor 
calculates  the  one-dimensional  normalized  MTF  of  the  focal  plane  array  as  a 
function  of  spatial  frequency  using  the  predicted  voltage  response  to  a  line 
source  input  as  described  in  Section  4.2.  This  calculation  involves  a  numerical 


integration  of  the  voltage  over  the  line  of  detectors  set  up  in  the  preprocessor 
for  this  calculation.  The  MTF  for  the  imperfect  charge  transfer  processes  in  the 
CCD  are  then  evaluated  and  their  product  formed  for  determining  the  total 
MTF. 

The  postprocessor  then  evaluates  the  sum  of  the  squares  of  the  various 
noise  currents  as  described  in  Section  4.2,  and  numerically  integrates  this  sum 
over  the  spatial  frequency  range  of  interest.  The  evaluation  of  the  noise 
currents  is  a  simple  matter  of  evaluating  algebraic  equations  given  in  Section 
4.2.  The  detector  responsivity  is  then  evaluated  from  the  zero  frequency  MTF 
divided  by  the  incident  radiant  power.  Finally,  the  NEFD,  NEP  and  D*  are 
calculated  from  the  sum  of  the  noise  and  the  responsivity  and  other  available 
detector  parameters. 


SECTION  6.  CODE  VALIDATION 


Validation  of  the  TFPA  code  proceeded  by  running  check  cases  whose 
results  were  either  compared  quantitatively  to  results  from  analytic  expressions 
or  studied  qualitatively  for  anticipated  results.  Check  Cases  A  through  D 
involved  the  heating  of  detector  slabs  by  either  a  constant  temperature  boundary 
or  a  constant  irradiance.  Check  Cases  E  through  J  had  detector  materials 
heated  by  modulated  irradiances.  All  check  cases  and  their  results  are  described 
below. 

CHECK  CASES  A  THROUGH  D 

Check  Cases  A  through  D  modeled  detector  slabs  whose  temperature 
changes  were  caused  by  constant  irradiances  or  constant  temperature  boundaries. 
The  temperatures  within  the  slabs  were  compared  to  temperatures  generated 
from  theoretically  produced  analytic  expressions.  Each  case  is  discussed  below. 

Check  Case  A 

The  initial  check  case  considered  the  constant  heating  of  a  single  detector 
whose  temperature  was  assumed  uniform  at  all  times  (the  so-called  "lumped 
approximation").  The  radiative  and  convective  heating  were  neglected  and  the 
pedestal  and  insulation  were  assumed  to  have  negligible  thermal  conductivity  so 
that  the  detector  was  effectively  insulated.  The  detector  was  heated  by  a 
constant  irradiation  E0  incident  on  the  top  surface. 

The  temperature  T  for  this  detector  as  a  function  of  time  (t)  was 
calculated  by  the  code  and  compared  to  the  theoretical  temperatures  given  by 
the  analytic  expression: 


T(t)  =  T(0)  +  E0At/CV 
where:  T(t)  =  Detector  temperature  (Kelvin) 


t  =  Time  (seconds) 

E0  =  Mean  irradiance  (w/m2) 


A  =  Irradiated  surface  area  (m2) 

C  =  Volumetric  heat  capacity  (J/m2-K) 

V  =  Detector  volume  (ma) 

The  temperature  rise  will  be  in  K  when  Eo  *s  *n  watts/m2.  This  expression  shows 
that  for  a  constant  irradiance  E0,  the  change  in  the  temperature  of  the  detector 
should  be  linear  with  time.  Run  A  showed  that  the  code  produced  the  same 
temperature  and  time  relationship  that  theory  predicts. 

Using  typical  dimensions  and  thermal  properties,  a  one-noded  detector  was 
heated  by  a  constant  0.05  w/m2  irradiance  for  1  msec  in  10  ps  time  steps.  The 
code  generated  the  identical  linear  temperature  versus  time  relationship  as  did 
the  theoretical  expression  (see  Figure  7). 

Check  Case  B 

Check  Case  B  introduced  one-dimensional  conduction  by  considering  the 
thermal  relaxation  of  the  detector  material  which  was  assumed  to  be  uniform 
initially  at  temperature  T0.  Then  the  temperature  of  the  bottom  surface  was 
suddenly  raised  to,  and  held  steady  at,  another  temperature  Tf.  Transverse 
conduction  was  neglected  by  using  perfect  insulation  around  the  detector.  The 
top  surface  was  also  perfectly  insulated  by  setting  the  incident  irradiation, 
radiative  losses,  and  convective  losses  to  zero.  The  exact  solution  for  the 
comparison  of  code  generated  temperatures  for  this  problem  is  given  by  Carslaw 
and  Jaeger  [11]  as: 


(78) 


where  K  is  the  thermal  conductivity  (w/m-K),  t  is  time  (seconds),  L  is  the  total 
thickness  (meters),  and  z  is  the  vertical  distance  (meters)  from  the  z=0  boundary 
at  the  top  to  the  point  of  interest  within  the  detector.  Run  B,  which  is 
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FIGURE  7.  TEMPERATURE  HISTORY  OF  AN  INSULATED,  THERMALLY  LUMPED  DETECTOR 
(HEAT  CAPACITY  OF  1.32x10-6  j/|Q  UNDER  CONSTANT  IRRADIANCE 
(0.05  W/M2)  ON  AREA  OF  1.03xl0"8  M2,  CHECK  CASE  A. 


described  below,  produced  temperatures  in  excellent  agreement  (see  Figure  8) 
with  those  predicted  by  the  exact  expression  above. 


In  Code  Run  B,  the  pedestal  and  detector  were  given  identical  dimensions 
and  thermal  properties  and  divided  into  ten  layers  each  so  that  a  detector  slab 
with  20  layers  was  modeled.  The  initially  isothermal  detector  was  set  upon  a 
constant,  slightly  elevated,  temperature  boundary,  the  substrate.  The  detector 
was  heated  from  the  constant  temperature  substrate  for  1  ms  in  10  us,  time 
steps.  The  code  generated  temperatures  that  are  in  excellent  agreement,  (i.e., 
indistinguishable  from  theory  on  scale  shown,  see  Figure  8)  with  those  predicted 
by  the  exact  expression  above. 
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COMPARISON  OF  THERMAL  FOCAL  PLANE  ARRAY  MODEL  WITH  EXACT  THEORY  (CASE  B). 


Check  Case  C 

Check  Case  C  again  considered  one-dimensional  transient  heat  conduction 
in  a  detector  slab,  but  for  this  case  the  heat  source  was  a  constant  irradianee 
on  the  top  surface  instead  of  a  constant  temperature  boundary  at  the  bottom 
surface.  In  this  case,  the  detector  bottom  and  sides  were  perfectly  insulated  and 
no  heat  loss  by  radiation  or  convection  was  allowed. 

Carslaw  and  Jaeger  also  give  the  exact  solution  for  this  problem  when  the 
initial  temperature  is  zero  Ill]: 


Very  good  comparisons  between  the  analytic  and  code  temperatures  resulted 
from  the  code  run  described  below. 

In  this  Run  C,  a  detector  of  typical  dimensions  (see  Table  1)  was  divided 
into  ten  horizontal  layers,  and  was  surrounded  by  materials  with  zero  thermal 
conductivity  to  obtain  perfect  insulation.  The  detector  was  again  heated  by  a 
0.05  w/m2  irradiation  for  1  ms  in  10  ys,  time  steps.  The  calculated  temperatures 
for  the  layers  agreed  well  with  the  theoretical  temperatures  (see  Figure  9). 

Check  Case  D 

Check  Case  D,  like  Check  Case  C,  considered  the  one-dimensional 
transient  heat  conduction  in  a  detector  slab  heated  by  a  constant  irradianee  on 
the  top  surface.  Again  the  sides  were  perfectly  insulated,  but  unlike  Check  Case 
C,  the  bottom  surface  was  a  constant  temperature  boundary  rather  than  an 
insulated  one. 

Carslaw  and  Jaeger  [11]  again  give  the  exact  solution  for  this  problem 
when  the  initial  detector  temperatures  and  the  constant  boundary  temperature 


are  zero: 


TABLE  1.  CHECK  CASES  A  THROUGH  J 


CONSTANT 

UNITOKH  IRRADIANCE 


E  z 

T(z,t)  “ 

W  -[<^>2(r)(S]  •  •1”E2-"GXD]| 

(80) 

Again,  very  good  agreement  between  analytic  and  calculated  temperatures 
resulted  from  the  code  run  described  below. 

In  Code  Run  D,  the  detector  and  pedestal  were  given  identical  dimensions 
and  thermal  properties  and  divided  into  ten  layers  each  so  that  a  detector  slab 
with  20  layers  was  modeled,  as  in  Case  B.  The  initially  isothermal  detector  was 
set  upon  a  boundary  whose  temperature  remained  constant  at  the  detector's 
initial  value.  The  detector  was  surrounded  by  perfect  insulation  and  heated  by 
a  constant  0.05  w/m2  irradiance.  The  resulting  layer  temperatures  agreed  well 
with  those  predicted  by  theory  (see  Figure  10). 

CHECK  CASES  E  THROUGH  J 

The  code  validation  process  continued  with  Check  Cases  E  through  J  where 
the  heating  of  detector  material  was  caused  by  a  modulated  irradiance.  Check 
Cases  E  through  G  compared  code  predicted  and  theoretical  temperatures, 
resulting  from  a  spatially  modulated  irradiance,  in  long,  thin  slabs.  Check  Case 
H  studied  an  approximately  circular  detector  array  with  an  impulse  irradiance  on 
its  center.  Check  Case  I  considered  a  linear  detector  array  irradiated  with  a 
temporally  chopped  line  source.  The  temperature  information  from  this  case 
was  used  to  determine  thermal  responsivities  and  crosstalk.  Check  Case  J 
studied  a  cruciform  detector  array  that  had  an  impulse  irradiance  on  its  center 
and  had  no  insulation  material  between  detectors.  A  more  detailed  description 
of  these  check  cases  follows. 
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Check  Case  E 


This  test  case  modeled  a  long,  thin,  perfectly  insulated  detector  heated  by 
an  irradiance  which  was  spatially  modulated  as  a  periodic  step  function. 
Emission  of  radiation  eventually  allowed  the  detector  to  approach  thermal 
equilibrium.  The  temperatures  along  the  detector  in  near  thermal  equilibrium  as 
predicted  by  the  code  were  compared  to  the  temperatures  derived  analytically. 
The  analytic  expression  for  temperature  difference  between  the  detectors  and 
the  ambient,  AT,  as  a  function  of  position,  x,  at  thermal  equilibrium  is  given  as 
(41 


AT(x)  -  A  ]1  +  £  B  sin  Wx  (81) 

L  m*0  "  S  J 

where:  A  =  E0  /  4oT03 

_ 1 _ 

.  kd(2nrt-l)2ir2 
4eoT 

o 

T0  =  Initial  temperature  and  temperature  to  which  heated 
detector  emits  radiation  (Kelvin) 

k  =  Thermal  conductivity  (joules/sec/K/m) 

d  =  Detector  thickness  (meters) 

S  =  Length  of  irradiation  step  (meters) 

This  expression  assumed  that  the  irradiance  is  off  in  the  region  where  -S<x<0 
and  on  where  0<x<S,  and  periodic  every  length  of  2S.  This  expression  showed 
that  the  temperature  at  equilibrium  was  periodic  every  2S,  and  that  no 
significant  temperature  gradients  exist  at  x  =  ±  S/2.  Thus  the  code  used  a 
detector  with  perfect  insulation  at  x  =  ±S/2  to  model  an  infinitely  long  detector 
under  a  periodic  step  irradiance.  Six  different  runs  were  made  for  Test  Case 
E  and  are  described  below.  Run  E-5  gave  temperatures  that  agreed  extremely 
well  with  the  analytic  temperatures. 


Bm  =  n(2nHT) 


The  code  used  a  2.5  mm  long  and  5  jxn  thick  detector  divided  into  20  x- 
direction  nodes  that  are  initially  at  ambient  temperatures.  Half  of  the  detector 
was  irradiated  with  0.05  w/m^  for  10,000  time  steps  of  45  ys  duration  each.  At 
the  end  of  the  run,  thermal  equilibrium  had  not  yet  been  attained;  more  time 
steps  were  needed.  The  node  temperatures  did  show,  however,  that  the 
theoretical  distribution  from  the  analytic  expression  had  begun  to  form  (see 
Figure  11). 

Run  E-2 

In  this  second  code  run,  the  analytic  temperatures  at  equilibrium  were 
inserted  as  the  code's  initial  temperatures  to  check  that  the  code  did  not  change 
them  during  the  run.  In  1,000  time  steps,  the  temperatures  had  been  changed 
negligibly  with  the  exception  of  the  temperatures  of  the  nodes  on  either  side  of 
the  on/off  boundary. 

Run  E-3 

From  Run  E-l  it  seemed  impractical  to  run  the  code  until  the  node 
temperatures  increased  from  ambient  to  equilibrium.  The  calculational  time 
step  of  45  ys  was  at  the  ceiling  of  allowable  increments  for  guaranteed  stable 
solutions,  (see  Section  4.1)  but  a  larger  time  increment  would  allow  for  a  faster 
approach  to  equilibrium  if  an  instability  would  not  cause  a  problem.  These  new 
runs  found  stable  solutions  for  4.5  ms  increments  but  unstable  for  increments  of 
1.0  seconds  or  greater. 

Run  E-4 

Node  temperatures  successfully  increased  from  ambient  to  their  respective 
equilibrium  temperatures  in  only  1,000  time  steps  when  45  ms  steps  were  used. 
These  equilibrium  temperatures  agreed  quite  well  with  the  analytic  tem¬ 
peratures,  except  at  nodes  at  the  on/off  boundary,  where  slight  differences 
occurred  (see  Figure  12). 
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FIGURE  11. 


CODE  RESULTS  FOR  CHECK  CASE  E,  SHOWING  THE  BEGINNING  OF  THE 
FORMATION  OF  THE  THEORETICAL  DISTRIBUTION,  RUN  E-l. 


RELATIVE  TEMPERATURE 


Run  E-5 


This  run  used  a  12.5  mm  long  and  5  ym  thick  detector  divided  into  60  x- 
direction  nodes  initially  at  ambient  temperature.  Again,  half  of  the  detector 
was  irradiated  with  0.05  w/m^  for  1,000  time  steps  of  45  ms  duration. 

The  equilibrium  temperatures  were  approached  more  quickly  in  this  run 
with  more  nodes,  and  the  temperatures  agreed  extremely  well  with  those  from 
the  analytic  expression  (see  Figure  13). 

Run  E-6 

The  previous  runs  used  a  detector  model  whose  only  means  of  losing  heat 
was  through  emission  of  radiation.  In  order  to  check  the  code’s  ability  to 
transfer  heat  by  convection,  an  equivalent  convection  coefficient  was  used  to 
replace  the  emissive  heat  loss  mechanism.  Thus  the  emissivity  was  set  to  zero 
and  the  convection  coefficient  was  chosen  to  give  identical  theoretical 
equilibrium  temperature  as  before,  using  the  same  detector  with  20  nodes. 

Again,  the  1,000  45.0  ms  time  steps  were  not  sufficient  to  allow  the  node 
temperatures  to  equilibrate.  Since  radiation  and  convection  are  dependent  on 
different  powers  of  temperature,  the  rate  at  which  the  two  heat  loss 
mechanisms  equilibrate  differ  (the  convective  power  of  one  being  slower  than 
the  radiative  power  of  four).  (Later  versions  of  the  code  use  convective  heating 
proportional  to  temperature  difference  raised  to  the  5/4  power.)  The 
temperatures  from  this  convective  loss  run  did,  however,  begin  to  show  the  same 
theoretical  temperature  distribution  as  the  radiative  loss  run  showed  (see  Figure 
14). 


Check  Case  F 

This  test  case  was  very  similar  to  Run  E  in  that  it  modeled  a  long,  thin, 
perfectly  insulated  detector  heated  by  a  spatiallv  modulated  irradiance.  In  this 
case,  however,  the  spatial  modulation  was  p  sinusoid  function  rather  than  a  step 
function.  This  case,  too,  allowed  heat  to  be  lost  by  the  detector  through 
emission,  so  that  thermal  equilibrium  eventually  resulted  after  heating  for  some 
time  period. 
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CODE  PREDICTED  NODE  TEMPERATURES 
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FIGURE  14. 


COMPARISON  OF  TEMPERATURES  FOR  HALF-IRRADIATED  SLAB  LOSING 
HEAT  BY  CONVECTION.  CHECK  CASE  E,  RUN  E-6,  AND  FROM  THEORY 
FOR  THE  EQUILIBRIUM  CASE.  EQUILIBRIUM  NOT  REACHED. 
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These  equilibrium  temperatures  along  the  detector  were  compared  to  the 
theoretically  predicted  temperatures  expressed  analytically  as: 

AT(x)  -  A  Jl  +  B  *in  [(2irnx  +  */2)]  j  (82) 


where: 


B  = 


1  + 


2 

n  w  ltd 
3 

eoT 


T0  =  Ambient  temperature  (Kelvin) 
n  =  Spatial  frequency  (cycles/meter) 
d  =  Detector  thickness  (meters) 


This  expression  shows  that  the  temperatures  are  periodic  every  l/n  distance  and 
that  no  appreciable  temperature  gradients  exist  at  x=0  and  x=l/n.  Thus,  the 
code  used  a  detector  with  perfect  insulation  at  x=0  and  x=l/n  to  model  an 
infinitely  long  detector  under  a  spatially  sinusoidal  irradiance. 

Four  runs  were  performed  for  Check  Case  F  as  described  below.  Run  F- 
3  gave  temperatures  at  equilibrium  that  agreed  very  well  with  the  analytic 
temperature. 

Run  F-l 

The  run  used  a  12.5  mm  long  by  5  vm  thick  detector  divided  into  20  nodes 
along  its  length.  Starting  at  ambient  temperature,  the  nodes  were  heated  under 
an  irradiance  which  varied  sinusoidally  along  the  detector's  length  from  zero  to 
0.1  w/m2  at  a  frequency  of  80  cycles/meter.  The  heating  occurred  for  10,000 
45  ys  time  steps. 

As  in  Run  E-l,  not  enough  time  was  used  for  the  nodes  to  reach  thermal 
equilibrium,  but  the  temperatures  were  tending  toward  the  analytically  expected 
distribution. 
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Run  F-2 


As  in  Run  E-2,  the  theoretically  predicted  node  temperatures  at  equili¬ 
brium  were  inserted  into  the  code's  initial  temperatures  to  determine  if  the  code 
would  change  them.  In  5,000  45  vs  time  steps,  no  significant  changes  were  made, 
implying  that  the  code  agreed  with  the  equilibrium  temperatures. 

Run  F-3 

It  was  found  in  Run  E-5  that  a  calculational  time  increment  which  could 
have  given  an  unstable  solution  actually  gave  stable  results.  So,  this  run  also 
used  1,000  45  ms  time  steps  to  attempt  the  heating  of  the  detector  from 
ambient  to  equilibrium  temperatures.  The  5  ym  thick  detector  was  divided  into 
60  nodes  along  its  12.5  mm  length  and  was  heated  by  80  cycles/meter  irradianee 
that  varied  from  0.0  w/m2  to  0.1  w/m2.  This  run  allowed  thermal  equilibrium 
to  be  reached  and  the  node  temperatures  agreed  very  well  with  those  predicted 
by  theory  (see  Figure  15). 

Run  F-4 

This  run  was  identical  to  the  previous  run  except  that  transverse 
conduction  was  allowed  through  a  common  electrode  instead  of  through  the 
detector.  The  product  of  the  thermal  conductivity  and  thickness  for  the 
electrode  in  this  run  was  identical  to  the  product  for  the  detector  in  the  last 
run,  so  that  identical  temperatures  for  the  two  runs  resulted. 

Check  Case  G 

Test  Cases  F  and  G  differ  only  by  the  frequency  of  their  spatially 
sinusoidal  irradiances.  Test  Case  G  modeled  a  thin,  perfectly  insulated, 
infinitely  long  detector  heated  by  an  irradianee  which  varied  sinusoidally  over 
the  detector's  length.  Upon  heating,  the  detector  radiation  increased  and 
eventually  the  detector  approached  thermal  equilibrium  where  the  rate  of 
emission  equaled  the  rate  of  absorption.  The  temperature  along  the  detector  at 
thermal  equilibrium  was  predicted  theoretically  by  the  relation  found  under  Case 
F,  which  serves  as  a  basis  of  comparison  for  this  code  check  case  also. 
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FIGURE  15.  COMPARISON  OF  EQUILIBRIUM  TEMPERATURES  FOR  DETECTOR  SLAB 
IRRADIATED  SINUSOIDALLY  ALONG  ITS  LENGTH  FROM  CHECK  CASE 
F,  RUN  F-3,  AND  FROM  THEORY. 


The  following  run  gave  equilibrium  temperatures  in  agreement  with  the 
analytically  produced  temperatures. 

Run  G-l 

A  5  p  thick  detector  was  divided  into  40  nodes  along  its  12.5  mm  length 
and  was  initially  set  at  ambient  temperature.  The  detector  was  heated  by  a 
spatially  sinusoidal  irradiance  which  varied  between  zero  and  0.1  w/m2  at  a 
frequency  of  800  cycles/meter.  The  heating  lasted  for  1,000  45  ms  time  steps. 

The  node  temperatures  in  this  run  successfully  approached  equilibrium  and 
were  in  agreement  with  those  temperatures  predicted  analytically. 

Check  Case  H 

Test  Case  H  modeled  a  small,  approximately  circular,  array  of  detectors 
receiving  a  short  impulse  of  irradiance  at  the  center  of  the  center  detector  (see 
Figure  16).  It  was  anticipated  that  the  array  would  transfer  the  heat  radially 
with  the  effective  temperatures  of  the  detectors  being  symmetric  about  any  axis 
that  passed  through  the  center  of  the  array.  This  run  showed  that  the  code 
handled  this  situation  properly. 

A  small,  roughly  circular,  detector  array  was  formed  from  37  detectors  and 
briefly  irradiated  at  the  center.  The  square  detectors  were  0.001  inches  on  a 
side,  divided  into  25  nodes,  and  set  on  square  pedestals  with  nine  nodes.  Typical 
thermal  properties  were  used  for  the  detectors,  pedestals,  and  insulations.  The 
irradiance  was  a  10  us  impulse  of  0.05  w/m2  on  the  center  node  of  the  array. 
The  code  gave  perfectly  symmetric  temperature  distributions  for  all  time  steps 
as  was  anticipated. 

Check  Case  I 

In  Test  Case  I,  one  detector  in  a  linear  detector  array  received  square 
wave  chopped  irradiation.  The  chopping  effect  caused  the  effective  temperature 
of  the  detectors  to  fluctuate  between  high  and  low  fixed  values  when  cyclic 
equilibrium  was  reached.  The  difference  between  these  high  and  low  tem¬ 
peratures,  the  effective  temperature  fluctuation,  gave  a  means  for  measuring 


FIGURE  16.  DETECTOR  ARRAY  MODELED  IN  TEST  CASE  H. 

the  thermal  responsivity  of  the  irradiated  detector  and  thermal  crosstalk 
between  the  irradiated  detector  and  any  other  detector. 

The  test  runs  for  this  case  were  first  used  to  validate  and  optimize  the 
method  of  modeling  a  linear  detector  array.  Later  runs  were  used  to  simulate 
the  effects  of  varying  detector  parameters  and  irradiances  on  the  effective 
temperature  fluctuations.  A  description  of  the  test  runs  follows. 

Run  1-1 


The  first  run  was  used  simply  to  verify  that  the  code's  linear  array 
responded  to  a  chopped  source.  Twenty  detectors,  divided  into  five  nodes  along 


their  lengths,  were  placed  into  the  array.  The  detectors  had  typical  TGS 
thermal  properties  (see  Appendix  A),  and  were  insulated  and  supported  on 
pedestals.  An  irradiance  of  0.05  w/m^  fell  on  the  first  detector  at  a  frequency 
of  about  500  Hz.  The  effective  detector  temperatures  that  resulted  did 
fluctuate  as  anticipated  and  the  first  few  detectors  reached  cyclic  equilibrium 
in  about  ten  chop  cycles. 

Run  1-2 

For  this  run  the  chopped  source  narrowed  to  only  one  node  width  and  the 
irradiated  detector  was  at  the  center  of  a  19  detector  linear  array,  and  the  chop 
frequency  was  reduced  to  around  50  Hz.  The  resulting  effective  temperatures 
were  symmetric  as  anticipated,  but  the  computational  time  was  excessive.  This 
run  required  approximately  13  hours  on  our  DEC  VAX  11/730. 

Run  1-3 

The  previous  run  gave  duplicate  effective  temperatures  on  either  side  of 
its  center  of  symmetry  and  thus  wasted  computational  time.  In  order  to 
consider  only  one-half  of  the  array,  however,  the  detectors  had  to  be  divided 
into  an  even  number  of  nodes,  with  the  two  center  nodes  being  irradiated  on  the 
center  detector.  This  gave  an  integral  number  of  nodes  on  either  side  of  the 
array's  center  of  symmetry.  The  run  here  described  served  as  a  check  for  a 
later  run  and  was  similar  to  the  last  run  in  all  respects  except  that  there  were 
six  detector  nodes,  two  nodes  being  irradiated  on  the  center  detector.  It  too 
gave  symmetric  effective  temperatures  around  its  center  of  symmetry  and  used 
excessive  calculational  time. 

Run  1-4 

This  series  of  runs  first  checked  the  code's  ability  to  work  with  only  one- 
half  of  a  symmetrically  irradiated  linear  array  and  then  studied  the  effect  of 
varying  irradiance  levels  on  the  effective  temperatures.  Only  ten  of  the  19 
detectors  of  the  previous  run  were  used,  but  all  other  array  properties  were 
unchanged  for  the  first  trial  run. 


The  ten  detectors  for  the  first  trial  run  behaved  identically  to  their 
counterparts  in  Run  1-3.  Thus  much  less  computational  time  was  required  for 
a  symmetric  case.  The  other  trial  runs  in  this  group  showed  that  the  effective 
temperature  fluctuation  at  cyclic  equilibrium  was  proportional  to  the  irradiance 
over  a  wide  range  of  values.  From  a  low  irradiance  value  of  0.025  w/m^  to  a 
value  of  75  w/m^  (corresponding  to  an  8-12  ym  radiation  from  a  target  whose 
temperature  exceeds  the  detector  temperature  by  about  100K),  the  effective 
temperature  variation  was  linear  with  irradiance. 

Run  1-5 

In  this  series  of  trial  runs,  linear  detector  arrays  consisting  of  different 
materials  were  irradiated  at  different  chopping  frequencies  in  order  to  study  the 
effects  on  the  thermal  responsivity  and  the  thermal  crosstalk  between  detectors. 
The  results  were  then  compared  to  similar  results  generated  by  a  HEOD  model 
[81.  A  detailed  description  of  results  from  these  runs  and  the  comparisons  is 
found  in  Section  7. 

Run  1-6 

This  series  of  runs  checked  to  see  how  the  width  of  the  "delta  function" 
of  irradiation  on  a  detector  in  a  linear  array  affected  the  thermal  responsivity 
and  the  thermal  crosstalk.  It  was  found  that  narrowing  the  irradiation  strip  from 
one-third  to  one-sixth  of  the  detector  width  resulted  in  negligible  changes  in 
detector  responsivity  and  crosstalk.  Widening  the  strip  to  include  all  of  the 
detector  surface,  however,  resulted  in  more  crosstalk  and  significantly  increased 
responsivities. 

Check  Case  J 

This  test  case  modeled  a  five  detector  cruciform  array  with  no  solid 
insulation  between  detectors.  This  array  was  given  a  short  irradiation  pulse  at 
the  center  node  to  determine  if  the  effective  temperatures  remained  symmetric 
about  the  center.  Unlike  previous  cases  where  heat  was  conducted  between 
detectors  through  insulation,  the  heat  transfer  between  detectors  in  this  case 
was  by  radiation.  Also  in  this  case,  heat  loss  by  convection  from  surfaces 
previously  contacted  by  insulation  was  included.  The  following  description  of 
trial  runs  describe  how  the  model  performed. 


This  trial  run  sent  a  20  ys  pulse  of  a  0.05  w/m^  irradiation  to  the  center 
of  a  cruciform  array.  The  array  consisted  of  five  25.4  ym  square  detectors 
separated  by  2.54  ym.  The  center  detector  was  heated  and  lost  heat  properly 
through  convection  and  conduction,  but  did  not  transfer  any  appreciable 
radiation  to  adjacent  detectors  (i.e.,  less  that  10E-38  w,  the  limit  of  the 
computer).  According  to  this  model,  the  detector's  temperature  must  be 
significantly  higher  before  appreciable  emissive  losses  occur. 


SECTION  7.  EVALUATION  OF  TYPICAL  THERMAL  FOCAL  PLANE  ARRAYS 


The  TFPA  code  was  used  to  model  several  linear  arrays  and  materials 
being  considered  by  the  community.  The  materials  chosen  for  study  consisted 
of  LiTaOg,  PVF2,  DTGFB  and  BaTiC>3,  which  are  representative  of  work 
currently  being  pursued  in  industry  [12].  A  quantitative  description  of  the  arrays 
is  given  below  followed  by  a  discussion  of  the  responsivities  and  crosstalks 
predicted. 

LiTa03,  PVF2  and  DTGFB  Detector  Arrays 

The  LiTa03,  PVF2  and  DTGFB  arrays  consisted  of  thin  sheets  of 
pyroelectric  materials  sectioned  into  50  pm  square  detectors.  The  detector 
thickness,  25  pm  for  LiTaC>3  and  DTGFB  and  6pm  for  PVF2,  was  uniform  over 
the  detector  area  except  for  a  one-half  thickness  notch  on  the  5  pm  border  of 
the  DTGFB  detectors.  Each  detector  contacted  the  electrode  beneath  it  by 
means  of  an  epoxy  pedestal  which  was  about  two-thirds  of  the  detector  area. 
The  LiTa03  array  was  surrounded  by  air,  while  the  other  arrays  were  evacuated. 

Each  linear  array,  built  with  a  central  detector  and  nine  more  on  either 
side,  was  irradiated  at  the  center  by  an  approximate  line  source  of  3.4x10"^ 
w/m2.  This  irradiance  represents  a  target  in  the  8-12  Pm  spectral  band,  which 
is  10K  warmer  than  ambient  (taken  to  be  298K).  The  irradiance  was  chopped 
at  50,  100,  200  and  400  Hz  in  order  to  study  the  thermal  responsivity  of  the 
center  detector  and  the  thermal  crosstalk  between  the  first,  second,  third  and 
fourth  detectors. 

The  code  was  usually  instructed  to  divide  the  individual  detectors  into 
eight  nodes  in  the  direction  of  the  array;  two  5  pm  nodes  at  the  borders  and  six 
nodes  for  the  remaining  40  Pm  length.  In  the  direction  perpendicular  to  the  line 
array,  only  three  nodes  were  used;  two  were  5  pm  wide  at  the  boundary  and  one 
was  40  pm  in  width.  For  each  frequency  studied,  the  code  operated  for  two  and 
one-half  chop  cycles,  except  for  the  400  Hz  case  which  ran  for  four  and  one- 
half  chop  cycles.  Irradiation  during  the  "on"  half  period  was  incident  on  only  the 
two  centermost  nodes  of  the  central  detector. 


BaTiC>3  Detector  Arrays 

The  BaTiOj  arrays  used  in  this  evaluation  consisted  of  distinct  25  ym 
square  detectors  that  were  12.7  y  m  thick.  Each  detector  was  positioned  on  an 
epoxy  pedestal  which  was  also  12.7  ym  thick,  but  had  only  two-thirds  of  the 
detector  area  in  contact.  The  detectors  were  separated  by  0.5  ym  and  a  glass¬ 
like  insulation  filled  the  gap  between  each  detector  pedestal  assembly. 

This  BaTiC>3  array,  built  with  a  central  detector  and  four  more  on  either 
side,  was  irradiated  at  the  center  by  an  approximate  line  source  of  2.16x10"'’ 
w/m2,  identical  to  the  irradianee  on  the  previous  three  arrays.  This  irradiance 
was  chopped  at  50,  100,  200  and  400  Hz  to  again  study  the  thermal  responsivities 
and  crosstalk. 

The  code  divided  the  detectors  into  nodes  in  a  slightly  different  manner  for 
the  BaTiOj  array.  In  the  direction  of  the  linear  array,  six  nodes  made  up  the 
detector  length  of  25  ym  and  a  2.5  yn  insulation  node  was  added  on  either  side. 
Perpendicular  to  the  array,  one  node  filled  the  25  ym  detector  width  and  a  5  ym 
insulation  node  was  added  on  either  side.  The  detector  and  epoxy  pedestal  each 
contained  one  node  vertically.  Five  detectors  were  used  by  the  code  to  model 
the  nine  detectors  in  the  symmetric  array.  Only  the  two  centermost  nodes  of 
the  modeled  central  detector  received  the  chopped  irradiation,  and  this  occurred 
for  only  2-1/2,  or  4-1/2  for  the  400  Hz  case,  chop  cycles. 

RESPONSMTY  AND  CROSSTALK 

The  thermal  responsivity  and  crosstalx  of  four  linear  detector  arrays  has 
been  determined  through  the  use  of  the  TFPA  code  at  four  frequencies  of 
chopped  irradiation.  These  results  are  shown  in  Table  2. 

The  thermal  responsivity  is  given  by 

*  -  g  (83) 

where  AT  (Kelvin)  is  the  effective  temperature  fluctuation  of  the  irradiated 
detector,  or  the  difference  in  volumetric  average  temperature  at  its  maximum 


TABLE  2.  NRC  SIMULATION  GENERATED  THERMAL  RESPONSIVITIES  AND  CROSS  TALKS 
FOR  TYPICAL  THERMAL  FOCAL  PLANE  ARRAYS. 


MATERIAL 

CHOPPING 

FREQUENCY 

(Hz) 

CYCLIC 

TEMPERATURE 

DIFFERENCE 

U) 

CROSS  TALK  (ACTUAL/ I APPARENT))* 

(*) 

(°K/W) 

DETECTOR  #2 

DETECTOR  13 

DETECTOR  14 

UTaOj 

50 

0.217E-4 

7,990 

54.8/(54.8] 

21.6/(21.61 

10.9/(6.31 

100 

0.157E-4 

5,780 

42.7/(42.7) 

15.7/19.6) 

4. 1/(0. 7) 

200 

0.106E-4 

3,900 

24.5/(  23.3) 

9. 1/(4. 0) 

5.2/13.51 

400 

0.700E-5 

2,580 

22.8/(16.61 

0.0/12. 7) 

0.0/13.4] 

pvf2 

50 

0.102E-3 

37,540 

3. 9/13.9) 

0.0/10.0) 

0.0/I0.0J 

100 

0.906E-4 

33,350 

2.0/10.7] 

o.o/to.o) 

0. 0/(0. 01 

200 

0.739E-4 

27,200 

0.4/10.2] 

o.o/to.o)  1 

0. 0/(0. 01 

400 

0.522E-4 

19,210 

0.0/10.2] 

o.o/to.o] 

0.0/10.0) 

DTGFB 

50 

0.587E-4 

21,610 

21.3/(21.3) 

3.0/11.0) 

0.5/10.31 

100 

0.428E-4 

15,750 

13. 6/(9. 8) 

0-8/10.6) 

0.0/10.1) 

200 

0.288E-4 

10,  600 

7.7/13.7] 

0.0/11.1] 

0.0/10.31 

400 

0.176E-4 

6,480 

3. 4/(1. 3) 

0.0/11.2) 

0.0/10.2) 

BaTIOj 

50 

0.185E-4 

16,890 

22.6/1  24.1] 

3.6/13.8] 

0.6/10.6] 

100 

0.173E-4 

15,790 

23.9/(  23.9) 

3.6/13.61 

0.6/10.51 

200 

0.162E-4 

14.790 

20.5/(20.5] 

2.5/1 1.6] 

0.3/10.11 

400 

0.130E-4 

11,870 

12.3/(11.3) 

1. 0/(0. 6] 

0.0/10.0) 

•IF  ACTUAL  IS  LESS  THAN  APPARENT;  CYCLIC  EQUILIBRIUM  WAS  NOT  APPARENT. 

NOTE:  IRRADIATED  AREA  -  0.533E-9  M2 

E  •  5.0941  W/H2  (CORRESPONDING  TO  THE  NET  IRRADIANCE  FROM  A  TARGET  WHOSE  TEMPERATURE 
IS  10  DEGREES  ABOVE  THE  DETECTOR  TEMPERATURE) 

and  minimum  over  a  chop  cycle,  E  (w/m2)  is  the  irradiance,  A  (m2)  is  the 
irradiated  area,  and  R  (K/w)  is  the  thermal  responsivity. 

Thermal  crosstalk  is  given  by 

AT 

CT  =  AT*  X  (100%)  (n=2,3,4)  (84) 

where  ATn  (Kelvin)  is  the  effective  temperature  fluctuation  for  detectors  not 
irradiated,  AT  (Kelvin)  is  again  the  effective  temperature  fluctuation  for  the 
irradiated  detector,  and  CT  (percent)  is  the  thermal  crosstalk  between  the 
irradiated  and  n*h  detectors.  Table  2  includes  an  actual  find  apparent  crosstalk 


term.  The  actual  crosstalk  is  determined  from  the  maximum  effective 
temperature  fluctuation  of  a  non-irradiated  detector.  Apparent  crosstalk, 
however,  is  determined  by  the  differences  in  effective  temperature  at  the  start 
and  finish  of  a  chop  half-cycle  which  is  not  necessarily  in  phase  with  the  times 
of  maximum  and  minimum  temperatures.  (Apparent  crosstalk,  then,  indicates 
the  signal  dispersion  that  a  chopper-synchronized  CCD  would  read.) 

The  thermal  responsivity  and  crosstalk  were  found  to  vary  as  the  array  and 
irradiance  varied.  The  following  discussion  reports  how  the  effects  of  chopping 
frequency,  array  materials  and  dimensions,  code  variables,  approximate  line 
source  widths,  and  irradiation  levels  changed  responsivity  and  crosstalk  values. 

Chopping  Frequency 

As  expected,  thermal  responsivity  drops  off  as  the  chopping  frequency 
increases.  The  calculated  responsivity  is  not  constant  between  50  and  100  Hz, 
so  an  even  lower  chopping  frequency  would  be  needed  to  observe  this  leveling 
off.  Thermal  crosstalk  decreases  as  the  chopping  frequency  increases.  Thus  a 
tradeoff  is  shown  for  responsivity  and  spatial  resolution. 

Array  Materials  and  Dimensions 

The  best  responsivities  came  from  the  PVF2  array,  the  remaining  arrays 
ranked  in  decreasing  order  of  their  responsivities  were  the  BaTiOj,  DTGFB  and 
LiTaOg  arrays.  A  very  rough  relationship  from  the  results  showed  that  the 
responsivities  at  lower  chopping  frequencies  were  nearly  proportional  to  the  heat 
capacities  and  thicknesses  of  the  detectors  and  inversely  proportional  to  the 
detector  thermal  conductivities. 

The  arrays  ranked  in  order  of  increasing  percent  thermal  crosstalk  are  the 
PVF2,  DTGFB,  BaTiC>3  and  LiTaC^.  Again,  a  very  rough  relationship  exists 
where  the  thermal  crosstalk  is  nearly  proportional  to  the  thermal  conductivity 
of  the  material  joining  adjacent  detectors. 


Code  Variables 

Neither  the  responsivity  nor  the  thermal  crosstalk  values  changed  sig¬ 
nificantly  when  two  vertical  nodes  were  used  in  each  of  the  pedestal  and 
detector  instead  of  one  node  each.  Similarly,  doubling  the  number  of  nodes  in 
the  linear  array  direction  had  little  effect  (i.e.,  less  than  0.1  percent). 

Approximate  Line  Source  Widths 

Both  the  responsivity  and  crosstalk  increased  when  the  irradiation  was 
incident  upon  the  entire  detector  surface  as  opposed  to  only  about  one-third  of 
the  surface.  However,  when  irradiation  was  incident  on  only  about  one-sixth  of 
the  surface  as  opposed  to  about  one-third,  there  was  no  significant  change  in 
either  responsivity  or  crosstalk. 

Irradianee  Levels 

No  significant  changes  in  responsivity  were  found  when  irradianee  levels 
changed.  The  effective  temperature  fluctuations  were  found  to  be  proportional 
to  irradianee  over  a  wide  range  of  values  (i.e.,  from  0.025  to  75.0  w/m2). 

HEOD's  MODELS 

The  LiTa03,  PVF2  and  DTGFB  detector  arrays  modeled  by  the  TFPA  code 
were  based  upon  similar  arrays  used  in  HEOD's  code  [8].  The  thermal 
responsivities  generated  by  both  codes  are  shown  in  Figure  17.  It  can  be  seen 
that  the  HEOD  code  yielded  lower  responsivities  in  each  case,  due  probably  to 
different  modeling  techniques.  Percent  crosstalk  values  for  each  code's 
treatment  of  the  LiTa03  array  is  shown  in  Table  3.  HEOD's  values  are  slightly 
higher  for  the  crosstalk  to  the  first  adjacent  detector,  but  smaller  for  the  second 
adjacent  detector  (32  percent  smaller)  and  significantly  smaller  for  the  third 
adjacent  detector  (factor  of  about  11  percent  smaller). 


76 


£  104[ 


1  2 


•  PVf2  (6  P«) 
e  DTGFB  (25  um) 
n  LiTaOj  (25  m) 

—  HONEYWELL  (8) 

--  NICHOLS  RESEARCH 


10  20  50  100  200  500  1000 

FREQUENCY  (Hz) 


FIGURE  17.  COMPARISON  OF  THERMAL  RESPONSIVITY  VERSUS  CHOPPING  FREQUENCY 
PLOTS  FROM  HEOD'S  NODAL  ANALYSIS  CODE  AND  FROM  NRC'S  TFPA 
CODE. 


TABLE  3.  COMPARISON  OF  THERMAL  CROSSTALK  RESULTS  FOR  HEOD  [8]  AND 
PRESENT  MODELS 


FIRST  ADJACENT  DETECTOR 


FREQUENCY 


HONEYWELL  MODEL 


MAXIMUM 


PRESENT  MODEL 


50 

57.5 

100 

41.3 

200 

28.1 

400 

24.1 

MAXIMUM  -  CROSS  TALK  EVALUATED  FROM  MAXIMUM 
TEMPERATURE  DIFFERENCE 

APPARENT  -  CROSS  TALK  EVALUATED  FROM  MAXIMUM 
TEMPERATURE  DIFFERENCE  WITH  READOUT 
SYNCHRONIZED  WITH  CHOPPER 


77 


SECTION  8,  RESULTS 


Some  of  the  more  important  results  of  the  study  are  summarized  in  this 
section.  The  accuracy  of  the  code  compared  to  exact  analytical  solutions  is 
shown  in  Table  4  for  some  of  the  validation  cases  considered. 


It  is  seen  that  the  code  has  adequate  accuracy  for  the  purpose  of  relative 
performance  evaluation  of  focal  plane  arrays.  The  check  cases  were  selected  to 
test  various  parts  of  the  code  under  conditions  comparable  to  laboratory  or  field 
tests. 


The  present  thermal  model  predicts  responsivities  that  are  about  a  factor 
of  two  larger  than  those  predicted  by  the  HEOD's  model  for  their  assumed 
configuration  and  materials.  The  exact  configuration  and  properties  used  by 
HEOD  were  not  given  in  available  documentation  [8].  The  crosstalk  predicted 
by  the  present  model  is  within  about  5  percent  of  that  predicted  by  the  HEOD 
model. 


TABLE  4.  PERCENT  ERROR  IN  CODE  VERSUS  THEORY  PREDICTION  AT  AT  VERSUS 
TIME. 


Mc  -  M, 


X/L 

NODE 

600  vs 

1000  vs 

TOP 

0.025 

6 

i 

-6.4 

-1.4 

MIDDLE 

0.475 

87 

-0.13 

-0.04 

TOP 

0.05 

6 

0.18 

0.08 

MIDDLE 

0.45 

42 

0.21 

0.12 

TOP 

0.025 

6 

0.13 

0.08 

MIDDLE 

0.475 

87 

-0.34 

-0.11 

SECTION  9.  CONCLUSIONS 


A  fully  operational  thermal/electrical  computer  model  has  been  developed 
for  evaluating  the  performance  of  pyroelectric  thermal  detectors.  This  model 
allows  evaluation  of  rather  general  design  configurations  with  adjustable 
dimensions  and  properties  (thermal  and  electrical),  under  arbitrary  field  or 
laboratory  conditions.  The  computer  code  has  been  thoroughly  validated  against 
a  series  of  exact  solutions  for  a  variety  of  theoretical  situations.  As  far  as  was 
practical,  the  model  was  based  on  first  principles  formulations  with  reasonable 
simplifying  conditions.  The  code  was  designed  for  a  high  degree  of  "user- 
friendliness,"  as  well  as  for  modularity  to  facilitate  code  understanding  and 
modification. 

The  model  was  checked  against  the  one-dimensional  thermal  model  results 
published  by  HEOD  and  found  to  give  similar  results.  Detector  responsivity 
predicted  by  the  HEOD  model  are  about  a  factor  or  two  smaller  than  those  of 
the  NRC  model.  This  difference  is  probably  due  to  differences  in  the  model  or 
input  parameters.  The  predicted  crosstalk  of  the  two  models  is  surprisingly 
close;  within  about  five  percent. 
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SECTION  10.  RECOMMENDATIONS 


This  preliminary  development  effort  provides  a  useful  design  tool,  and 
uncovers  the  important  phenomenology  for  future  developments.  NRC  recom¬ 
mends  that  the  code  be  exercised  in  a  well-planned  series  of  computer 
"experiments"  to  evaluate  the  effects  of  various  design  parameters  on  the 
theoretical  performance  of  typical  sensor  concepts.  Several  code  upgrades  and 
modifications  are  desirable  to  investigate  more  accurately  the  key  phenomeno¬ 
logical  effects.  Specifically,  NRC  recommends  that  the  following  technical 
analysis  and  development  tasks  be  completed  as  soon  as  possible. 


Task  1.  Define  the  model  input  parameters  to  characterize  several  leading 
pyroelectric  detector  systems  of  interest  (e.g.,  systems  by  TI,  HEOD). 

Task  2.  Design  a  series  of  computer  experiments  to  evaluate  performance 
sensitivity  to  design  parameters  (e.g.,  materials,  dimensions,  con¬ 
figuration,  chop  rate,  readout  rate,  processing  noise  parameters). 
The  major  output  of  this  task  should  be  a  matrix  of  design  parameter 
values  over  which  the  sensitivity  analysis  would  be  run,  and  a  set  of 
nominal  values. 

Task  3.  Simulate  the  performance  of  the  matrix  of  sensors  (defined  in  Task 
2)  for  an  actual  irradiance  scene,  measured  by  a  high  performance 
calibrated  infrared  imaging  system  filtered  in  two  or  more  spectral 
bands  (e.g.,  3-5  ym  and  8-12  ym).  The  performance  sensitivity  should 
be  characterized  by  parameters  such  as  MTF  (thermal  and  electrical); 
thermal  crosstalk  to  adjacent  detectors  and  detectors  two  and  three 
away  from  an  irradiated  detector;  detector  responsivity  as  a  function 
of  electrical  frequency  (w);  NEP  and  NEAT  versus  w;  D*. 

Task  4.  Various  implicit  techniques  (e.g.,  Crank-Nicholson)  should  be  examined 
for  possible  application  in  the  TFPA  code.  If  reasonable,  the  code 
should  be  modified  to  use  an  implicit  technique  as  an  option,  and 
comparisons  made  with  the  existing  explicit  technique. 

Task  5.  The  TFPA  code  should  be  modified  to  consider  the  effects  of  heat 
conduction  through  the  bottom  electrode. 

Task  6.  A  more  accurate  treatment  of  the  radiative  heat  transfer  between 
adjacent  detectors  should  be  developed. 

Task  7.  Formulate  necessary  equations  and  expand  code  to  provide  the  system 
throughput  in  response  to  the  application  of  arbitrary  stationary 
radiation  patterns  to  the  focal  plane. 


Task  8.  Formulate  necessary  equations  and  expand  code  to  provide  the  system 
throughput  in  response  to  the  application  of  arbitrary  time-varying 
radiation  patterns  to  the  focal  plane. 

Task  9.  Investigate  system  MTF  as  a  function  of  geometrical  point  of 
application  of  6 -function  radiation  input;  perhaps  devise  more 
representative  way  in  which  to  formulate  system  MTF. 

Task  10.  Using  representative  system  and  operational  parameters,  investigate 
the  limits  of  strict  applicability  of  the  code  with  respect  to  linearity, 
appropriateness  of  the  use  of  MTF,  etc.  Investigate  correlation  of 
behavior  in  orthogonal  directions  on  focal  plane.  A  consideration  of 
alternative  system  figures-of-merit  which  are  appropriate  to  thermal 
focal  planes  might  be  useful. 

Task  11.  Develop  the  necessary  software  to  produce  computer  generated 
graphical  output  of  various  simulation  output  parameters  of  interest. 
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APPENDIX  A.  THERMAL  AND  ELECTRIC  PROPERTIES  OF  PYROELECTRIC  MATERIALS  AT  AMBIENT  TEMPERATURE.  Numbers 
in  brackets  refer  to  literature  survey  bibliography  (Reference  4). 


APPENDIX  B.  THERMAL  PROPERTIES  OF  PEDESTAL,  ELECTRODE  AND  SUBSTRATE  MATERIALS.  Numbers  in  brackets 
refer  to  literature  survey  bibliography  (Reference  4)  and  numbers  in  parentheses  refer 
to  references  in  this  report. 


APPENDIX  C 


In  many  radiation  detection  schemes  utilizing  a  two-dimensional  array  of 
detector  elements,  readout  of  the  array  is  accomplished  using  standard  CCD 
devices.  Discussions  of  the  performance  properties  of  these  devices  under  a 
wide  variety  of  incident  image  characteristics  have  been  published  by  several 
authors  (see,  for  example.  References  8  and  10).  For  the  operation  of  a 
pyroelectric  detector  array,  the  incident  radiation  must  be  modulated  because 
the  pyroelectric  charge  is  generated  by  polarization  changes,  which,  for  small 
temperature  differences,  are  proportional  to  the  change  in  temperature  of  the 
pyroelectric  material.  In  the  following  analysis,  the  Fourier  transform  of  the 
voltage  output  of  a  finite  pyroelectric  array  with  CCD  readout  is  derived.  In 
addition  to  utilizing  a  detailed  model  of  the  CCD  sampling  process,  the  action 
of  a  chopper  is  explicitly  treated.  The  general  result  for  the  Fourier 
transformed  array  voltage  output  as  a  function  of  electrical  frequency  is  quite 
complicated.  Hence,  to  illustrate  the  main  features  of  the  general  result  which 
is  given  in  Equation  28,  discussion  of  a  simple,  chopped,  line  source  is  included 
at  the  end. 


The  conversion  by  a  pyroelectric  detector  array  of  a  chopped,  static 
(constant  in  time)  irradiance  input  to  a  voltage  output  is  assumed  to  oe  a  linear, 
convolutionary  process.  As  such,  the  output  from  say  the  (m,n)th  detector  is 
representable  as 


,v»  ■  lil 


dxd^dn  H(xm-£;,yn-n,t-T)  i(^,n)  c(xmH-vxT>yn-ri-vyT) 


where:  H(x,y,t)  =  Pyroelectric  array  transfer  function  (including  transfer 

inefficiencies  of  the  readout  process) 

I(x,y)  =  Input  irradiance  on  array  surface 

C(x,y)  =  Chopper  geometry 

(xm,yn)  =  Focal  plane  coordinates  for  (m.n)**1  detector  (see  Figure 
1) 

(Vx,Vv)  =  Focal  plane  coordinate  components  of  the  chopper  velo¬ 
city. 


85 
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FIGURE  18.  FOCAL  PLANE  COORDINATES  FOR  (M,N)th  DETECTOR. 

Equation  85  is  applicable,  under  the  stated  conditions,  only  if  the  chopper 
contributes  a  negligible  background  compared  to  the  image  scene's  irradiance.  If 
this  is  not  the  case,  the  chopper  contribution  to  the  irradiance  on  each  detector 
can  be  accounted  for  by  writing 

K£,n)  c(x-5,y-n)  -  i(£,n)  c(x-c,y-n)  +  i  [l  -  c(x-f,y-n)]  (86) 

where  Ic  represents  the  constant  irradiance  from  the  chopper.  Note,  only  when 
the  chopper  is  "closed,"  i.e.,  C=0,  does  I0  irradiate  the  array  elements  beneath 
the  chopper.  When  the  right  hand  side  of  Expression  86  is  inserted  into  Equation 
85  above,  l  he  term  containing  Ic  alone  contributes  only  to  an  unimportant  dc  or 
zero  frequency  component  of  the  output  power,  hence  may  be  dropped.  The  two 
remaining  terms  can  be  combined  into  an  expression  identical  to  Equation  85 
with  I(o$)  referenced  to  Ic.  Hence,  without  loss  of  generality,  Equation  85  will 
be  taken  to  represent  the  voltage  output  of  the  (m.n)1^1  detector  including  the 
chopper  irradiance. 

The  action  of  the  chopper  in  this  model  is  assumed  to  be  uniform 
translation  in  the  positive  x-direction,  hence  Vy=0.  In  addition,  the  chopper 
width  is  assumed  to  be  an  integral  number  of  detector  lengths,  i.e.,  jAx,  j  = 
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integer  >0,  and  its  length  is  essentially  infinite  (in  the  y-direction).  A 
representation  of  this  chopper  is 


e(x  -C-v  t»y -n) 

m  x  m 


(87) 


where 


(88) 


At  t=0,  the  first  j  rows  of  detectors  are  covered  by  this  chopper. 

A  somewhat  simplified  readout  process  is  included  in  the  model.  Readout 
is  represented  as  a  sampling,  in  time,  of  the  voltage  output  of  Equation  85  with 
samples  being  read  out  sequentially  in  the  order 

OUTPUT 

♦ 

DETECTOR  (0.1) 
t 

DETECTOR  (1,1) 

♦ 

OETECTOR  (2,1) 

♦ 

♦ 

DETECTOR  (H-1,1) 
t 

OETECTOR  (0,2) 

T 

OETECTOR  (1,2) 

♦ 

OETECTOR  (2,2) 
t 

♦ 

OETECTOR  (A- 1.2) 
t 

ETC. 


for  a  total  of  MxN  discrete  spatial  samples  of  the  entire  array  in  one  frame 


To  establish  notation,  and  to  be  more  precise  about  the  spatial  and 
temporal  sampling,  write  the  output  from  one  (arbitrary)  point  (x,y)  in  the  focal 
plane,  taken  at  time  t,  as 

V(x,y,t)  0  <.  x  <  (M-l)Ax 

V(x,y;t)  -  1  <  y  <  NAy  (89) 

0  Otherwise 

V  given  here  by  Equation  85  with  (x,y)  replacing  (xm»ym^  Upon  sampling, 
spatially,  the  voltage  output  will  become  a  doubly  periodic  sequence,  i.e., 
periodic  in  x  with  period  M  and  periodic  in  y  with  period  N.  This  spatial  output 
can  be  represented  as  one  period  in  M  and  one  period  in  N  of  the  sequence 

OO  00 

V(x,y;t)  -  l  l  V(x+kM,y+£N;t).  (90) 

k*-» 

where  V  is  given  by  Equation  86.  This  V  has  the  proper  spatial  behavior  and 
represents  the  focal  plane  output  which  is  sampled  spatially  and  temporally  (see, 
for  example,  Reference  14). 

Spatial  sampling  is  accomplished  by  taking  2J  spatial  samples  over  one 
detector  length.  Ax,  at  one  detector  sample  per  frame  time.  The  sample  time 
interval  is  derived  from  the  sampling  theorem  for  band-limited  functions.  That 
is,  if  all  temporally  varying  inputs  are  band-limited  to  a  band  of  width  2/T  sec-*, 
optimal  sampling  is  at  T/2  second  time  intervals.  With  a  scan  velocity  Vx 
over  the  smallest  length  scale,  Ax/J,  T  =  Ax/JVx.  Spatial  sampling  in  the 
direction  of  chopper  motion  is  handled  similarly.  Without  loss  of  generality,  we 
assume  spatial  frequencies  no  greater  than  J/Ax  are  present  in  the  scene,  hence 
sampling  is  done  in  intervals  of  length  A  x/2J.  Sampling  in  y  results  simply  from 
indexing  the  temporal  samples  so  as  to  obtain  correspondence  between  the 
readout  times  and  the  detectors'  y-coordinates  (see  below). 


At  this  stage  of  our  detector  array  modeling,  steady-state  response  has 
been  the  focus  of  our  analysis.  Hence,  in  Equation  85,  the  system's  temporal 
response  will  be  assumed  from  here  on  to  be  instantaneous.  In  this  case,  the  dt 
integration  becomes  trivial,  reducing  Equation  85  to  the  form 

V(x,y;t)  =  /  /  dxldyl  Hfx-Xj.y-yj)  Cfx-Xj-V^t)  Kxj.yj) 

"  (91) 

=  H(x,y)  C(x-Vxt)  **  1 


where  **  is  used  to  signify  a  two-dimensional  convolution.  At  this  time,  we  will 
apply  the  sampling  function  which  yields  the  CCD  readout  protocol  described 
above.  This  sampling  function  incorporates  the  two  characteristic  clock  cycles 
inherent  in  our  readout  procedure.  The  fastest  of  these  clocks  is  the  CCD 
device  clock  which  runs  at  a  rate  such  that  an  entire  array  is  read  out  in  a 
frame  time,  ts  =  £k/2JVx.  For  our  MxN  array,  this  clock  has  a  period  ofts/MN 
seconds  and  in  the  sampling  function,  the  summation  index  j  is  associated  with 
this  clock  period.  Thus,  to  accomplish  array  readout  in  a  frame  time,  a  single 
detector  voltage  is  output  each  CCD  clock  cycle  beginning  with  the  (x  =  0  •  Ax, 
y  =  iy)  spatial  location,  i.e.,  the  upper  left  hand  detector  in  Figure  1  alone. 
During  the  next  CCD  clock  cycle,  detector  (x  =  Cx,  y  =  Ay)  is  readout,  etc.  After 
M  CCD  clock  cycles,  the  y-coordinate  is  incremented  and  the  x-eoordinate  is 
reset.  After  MN  CCD  clock  ticks,  one  frame  time,  t  is  incremented  one  frame 
time,  and  the  chopper's  spatial  location  is  incremented  by  Ax/2J  (the  chopper 
blade  moves  in  the  x-direction  only).  The  x  and  y  coordinates  are  reset  and  the 
readout  procedure  begins  for  the  next  frame.  With  V  of  Equation  91  obeying 
Equation  89  above,  the  spatially  and  temporally  sampled,  spatially  periodic  array 
output  voltage  at  time  t  can  be  expressed  in  the  form 


V(x,y,t)  -  l  H(x,y)  C(x-V  t)  **  I  5  |  x  - 

J— »  / 


(mn) 

- 2J -  +  1  B»od(M)  Ax 


•  5  [y  -  (£)  »od(N)Ay]  •  6  (t-T.fgJ]  ) 
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where  (a)'  means  take  the  integer  part  of  a  and  <5{  }is  the  usual  delta  function. 
Rewriting  Equation  92  in  a  more  convenient  form  by  regrouping  terms,  the 
integer  part"  notation  can  be  incorporated  to  yield 


v(x,y,t) 


»  M-l  N-l 

ll  l  H(x,y)  C(x-V  t)  **  1 

k-co  y— (M-l)  A—  (N-l) 


•  6 


x  -  nod^J)  +  (j+iHfkNM)  mod(M)j  Ax 


•  6  {y  -  (A+kN)  mod(N)&y}  6{t+kTg)  (93) 


where  the  general  result 


®  M-l  N-l 

l  l  I  A(  j+MH-kNM ,  A+kN ,  k) 

k»-<»  j— (M-l)  A—  <j»-l) 


(94) 


has  been  used.  The  sequence  of  terms  in  Equation  93  whose  indices  are 
expressed  modulo  an  integer  can  be  treated  as  a  single  period  of  an  infinite 
periodic  sequence.  For  example,  consider  the  k  Mod  (2J)  index  in  the  first  delta 
function  in  Equation  93.  Representing  a  general  term  in  the  summation  over  k 
as  Z  (k  Mod  (2J)]  where,  without  loss  of  generality,  all  other  k  dependence  may 
be  temporarily  neglected.  For  0  <.k  <.2J-1,  the  index  k  Mod  2J  is  simply  k, 
and  is  zero  for  k  outside  this  interval.  Thus  in  this  index,  k  Mod  (2J),  Z  (k 
Mod(2J)]  is  a  sequence  of  length  2J.  This  sequence  can  be  considered, 
alternatively,  as  a  single  period,  of  duration  2J,  of  an  infinite  periodic  sequence 

■r. 

Z  [k]  where 


00 

zlk]  -  l  Z  [k+2rJ] 


(95) 


Applying  this  to  the  indices  k  Mod(2J),  (j+tM+kMN)  Mod(M)  and  U+kN)  Mod(N)  of 
Equation  93,  the  following  expression  results  where  summation  indices  have  been 
combined  when  possible. 


®  M-l  N-l  •  • 

V(x,y,t)«  II  III  H(x,y)  C(x-Vxt)  **  I 

k«- °°  j— (M-l)  £»-(N-l)  ot»-“ 

•  6  |  x  -  +  a  +  pi)  Ax  | 

•  6  {y  -  (JH-kM)Ay) 

•  6{t-[j+fcM4>!Om]Ts}  (96) 

Prior  to  Fourier  transforming  in  space  and  time,  introduce  the  following 
expression  for  the  three  delta  distributions  in  the  last  expression: 


6(x) 


(97) 


Interchanging  the  order  of  integration  and  summation,  Equation  96  then  takes 
the  form: 


V(x,y,t) 
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The  summations  can  be  performed  using  the  result 


J  eing  «  2n  J  6(g-2im) 
n“-°°  n“-®° 


to  yield 
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Fourier  transforming  this  last  expression  in  (x,y,t)  yields,  upon  interchanging  the 
order  of  integration, 


oo 

-  M-l  N-l  fff 

V(k  ,k  ,w)  *  2  n  l  l  dS  c 

y  j  — (M-l)  fi. —  (N-l>  JJJ 


ds2ds3  I 


where: 
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y  I(x,y)  (104) 


where  H  is  the  array  system  transfer  function  including  the  geometry  (i.e., 
spatial  extent  of  array,  detector  length,  width),  C  is  the  chopper  and  1  is  the 
input  irradiance.  The  Si's  (spatial  frequencies)  are  defined  as 
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kx  and  ky  are  the  spatial  frequencies  of  the  input  irradiance  and  w  is  the  output 
electrical  frequency.  Performing  integrations  in  Equation  101, 
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The  summations  over  j  mid  A  can  be  done  to  yield 


for  the  Fourier  transform  of  the  array  output  voltage  as  a  function  of  spatial 
and  electrical  frequencies.  Note  that  ftx,  fly  and  Az  are  functions  of  a, 3 
and  6  (see  Equations  109  through  111).  This  expression  contains  the  general 
relationship  of  input  spatial  frequencies,  kx  and  ky,  to  the  output  electrical 
frequency,  w.  This  can  be  brought  out  by  inserting  the  Fourier  transform  of  the 
chopping  function,  the  system  transfer  function  and  an  irradiance  distribution. 
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Using  Equation  88  above,  the  transform,  C(^),  is: 
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where 


(115) 


For  a  system  transfer  function,  assume  a  simple  MxN  rectangular  array  (as  in 
Figure  85  above).  The  Fourier  transform  of  the  system  function  is  (up  to  a 
phase) 
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where  <KM,N)  is  an  inessential  phase  dependent  only  on  M  and  N.  To  specify  an 
irradiance,  I(x,y)  is  taken  to  be  a  line  source  of  infinite  extent  in  the  y-direction 
lying  across  the  pth  row  of  detectors  at  an  x-coordinate  of  (p  +  q/2J)  Ax,  i.e.. 


X(x,y)  -  (p  +fj)Axj  , 


p  *  0,1,. . . ,M-1 
q  “  1,2,. . . ,2J 


hence  its  Fourier  transform  is 
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I(x,y)  has  been  normalized  such  that  summing  over  a  row  of  detectors  in  the  y- 
direction  yields  the  total  irradiance  input  to  the  array.  Instead  of  using  these 
transforms  in  the  general  expression  (Equation  113),  Equation  109,  written  in  a 
slightly  different  form,  will  be  used.  To  simplify,  Equation  109  for  the  example 
being  considered,  using  the  delta  function  in  f2y  which  appears  in  the  example 
irradiance  (see  Equation  118  above),  the  sum  over  6  can  be  given  the  form 
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where  the  summation  index,  6,  has  been  replaced  by  Y.  In  addition,  the  Fourier 
transform  of  the  system  function  H(ftx,fy)  can  be  simplified  to 
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by  utilizing  the  fty=0  condition  from  Equation  118.  Further  simplification  can  be 
obtained  by  noting  the  summation  index  a  appears  only  in  the  combination  a- 
8M  with  the  index  3  (see  Equations  110  and  112).  Hence  shifting  the  summation 
in  a,  Equations  110  through  112  become 
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The  summation  in  8 

also  can  be 

shifted,  i.e.. 
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With  this  change, 


Using  these  simplifications,  Equation  107  is  now  of  the  form 
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Since  only  the  bracketed  terms  in  this  last  equation  are  dependent  upon  a,  the 
summation  over  a  can  be  performed.  Focusing  on  the  bracketed  terms  only,  the 
summation  in  a  is  of  the  form 
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where:  A 
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Using  Poisson's  summation  formula,  the  sum  in  Equation  128  can  be  carried  out 
to  yield 
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Using  the  restriction  in  Equation  130,  Equation  129  becomes 
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where  a  satisfies 


With  these  results,  Equation  128  becomes 


Before  using  this  in  Equation  127,  it  is  convenient  to  note  an  alternative  form 
for  C(f^)  appearing  in  that  equation.  6(^),  given  in  Equation  114  above,  can 

be  written  as 
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Using  this  form  for  &(^z)  and  the  results  from  Equation  136,  Equation  127 
becomes 
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Using  the  limits  04  from  Equations  134  and  135,  the  summation  in  a  can  be 
carried  out  yielding 
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where 


0  =  f  +  fj  <k,  -  kcX> 


and  A  is  defined  in  Equation  128.  For  M  large,  the  limit  as  M*»  of  Equation 
139  provides  a  reasonable  approximation  to  sampled  array  behavior.  Taking  this 
limit,  the  a  summation  becomes 
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Equation  138  now  becomes 
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For  consistency,  the  large  M,N  limit  should  be  taken  in  this  last  expression  as 
done  in  the  summation  over  a.  In  doing  this,  the  sums  over  j  and  &  diverge  as 
2N  and  2M,  respectively.  However,  ws  behaver  as  (MN)-1,  thus  together 
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Some  care  is  required  in  taking  the  large  N  limit  in  the  summation  over  the  ky 
delta  function.  This  limit  is 
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This  result  follows  by  first  applying  the  sampling  distribution  of  Equation  144  to 
test  functions  (x)  which,  without  loss  of  generality,  are  taken  to  be  normalized, 
i.e.. 
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Replacing  the  6  distributions  with  the  representation. 
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the  sampling  represented  by  Equation  144  takes  the  form 
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where  k0y  has  been  pulled  out  of  the  6  distribution  using  the  relation 
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so  x  =  ky/koy  is  dimensionless.  Changing  variables.  Equation  147  becomes 
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where  the  integral  has  been  written  as  a  sum  of  three  terms.  Since  <t>  (x  + 
y/N)x  1  js  integrable  over  the  (-“,  -A)  and  (A,®)  intervals,  it  follows  from  the 
Riemann-Lebesgue  lemma  that  the  first  two  terms  of  Equation  149  tend  to  zero 
as  t-**.  For  A  sufficiently  small,  <Mx  +  y/N)  ~  <t>  (Y/N)  since  4>  is  continuous.  With 
this  Equation  149  becomes 
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With  a  change  of  variables,  xt=n,  and  the  result 
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Equation  150  takes  the  form 
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For  large  N,  the  summation  can  be  replaced  by  a  Riemann  integral  with  the 
identifications 
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to  yield 
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which  is  independent  of  N,  hence  the  limit  becomes  trivial.  Since  the  test 
functions  were  taken  to  be  normalized,  the  result  on  the  right  hand  side  of 
Equation  144  follows.  Combining  this  and  the  result  in  Equation  143,  the  voltage 
output  of  Equation  142  becomes 
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The  summation  over  A  in  this  result  is  simply  the  Fourier  transform  (in  the 
frequency  variable  kx+kox9  of  the  square  pulse  sequence  depicted  in  Figure  19a. 
The  exponential  is  a  phase  shift  reflecting  the  chopper  initial  position  and  it  has 
no  effect  on  the  amplitude  spectrum.  These  pulses  are  of  width  j&,  i.e.,  the 
physical  chopper  blade  covers  J  detector  elements  each  of  length  Ax  in  the  x- 
direction,  the  direction  of  chopper  motion. 


The  spacing  between  pulses  is  2jAc,  so  the  chopper  has,  in  this  model,  a  50 
percent  duty  cycle.  At  t=0,  the  edge  of  a  blade  is  aligned  with  the  line  source 
of  irradiance  at  the  x  coordinate  (p  +  q/2J)&.  The  chopper  and  line  source  are, 
as  mentioned  earlier,  assumed  to  be  of  infinite  extent  in  the  y-direction  (out  of 
the  page  in  Figure  19).  The  Fourier  transform  of  this  pulse  train,  which  appears 
in  Equation  155,  is  shown  in  Figure  19b  for  3  equal  to  a  fixed  integer.  The 
transform  is  a  line  spectrum  [entered  at  multiples  of  k^  =  2  it/ ax,  each  line  of 
which  has  side  bands  of  decreasing  amplitude  at  multiples  of  1^=2  ir/j Ax)). 

In  Figure  19c,  the  voltage  amplitude  spectrum  is  plotted,  at  fixed  B,  as  a 
function  of  normalized  electrical  frequency  w/Vx  whose  units  are  radians  per 
length.  The  complete  voltage  amplitude  spectrum,  for  the  infinite  array- 
chopped  infinite  line  source  approximation,  would  be  a  superposition  of  an 
infinite  set  of  spectra  like  Figure  19c,  each  with  a  different  integer  value  B . 
This  superposed  set  of  spectra,  one  for  each  integer  value  ofS,  is  the  result  of 
the  specialization  of  the  general  equation.  Equation  113,  to  the  case  of  a 
sampled,  infinite  two-dimensional  array  with  a  chopped  line  irradiance  source. 

In  the  case  where  the  two-dimensional  array  is  treated  as  finite,  each  of 
the  pulses  appearing  in  Figure  19c  would  appear  as  broad  I  sin  Ax/xl  pulses.  This 
would  in  turn  potentially  cause  serious  aliasing  difficulties.  In  an  approach  to 
minimizing  such  problems  by  a  study  of  tradeoffs  involving  chopper  char¬ 
acteristics  as  well  as  detector  response  characteristics.  Equation  113  should 
prove  to  be  extremely  useful. 

When  fully  integrated  with  the  code  developed  during  this  contract,  the 
general  treatment  given  in  this  appendix  (expressed  in  Equation  113)  of  a 
sampled  two-dimensional  rectangular  array,  with  an  arbitrary,  chopped  ir¬ 
radiance  input,  should  provide  a  significant  tool  in  simulation  and  design  of 
pyroelectric  arrays. 
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(LOCATION  OF  LINE  IRRADIANCE  SOURCE) 

(a)  SCHEMATIC  OF  THE  SPATIAL  EXTENT  (IN  THE  DIRECTION  OF  MOTION)  OF 
THE  CHOPPER;  EACH  BLADE  IS  j  DETECTOR  ELEMENTS  WIDE,  AND  SPACED 
21  ELEMENTS  APART.  INITIAL  POSITION  OF  CHOPPER  BLADE  RELATIVE 
TO  THE  LINE  SOURCE  IS  ALSO  INDICATED. 


AMPLITUDE  DECREASES  AS  ± 
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(b)  SPECTRUM  OF  SAMPLED  SPATIAL  FREQUENCIES  FOR  3  EQUAL  TO  FIXED  INTEGER. 
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(c)  VOLTAGE  AMPLITUDE  SPECTRUM  AS  A  FUNCTION  OF  NORMALIZED  ELECTRICAL 
FREQUENCIES  (NORMALIZED  TO  CHOPPER  VELOCITY)  FOR  THE  SAME  VALUE 
OF  3  AS  APPEARS  IN  THE  SPATIAL  FREQUENCY  SPECTRUM  OF  (b)  ABOVE. 


FIGURE  19. 


